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Abstract 

Functionally graded materials (FGMs) are characterized by spatially variable microstruc- 
tures which are introduced to satisfy given performance requirements. The microstructural 
gradation gives rise to continuously or discretely changing material properties which compli- 
cate FGM analysis. Various techniques have been developed during the past several decades 
for analyzing traditional composites and many of these have been adapted for the analysis of 
FGMs. Most of the available techniques use the so-called uncoupled approach in order to an- 
alyze graded structures. These techniques ignore the effect of microstructural gradation by 
employing specific spatial material property variations that are either assumed or obtained by 
local homogenization. The higher-order theory for functionally graded materials (HOTFGM) is 
a coupled approach developed by Aboudi et al. (1999) which takes the effect of microstructural 
gradation into consideration and does not ignore the local-global interaction of the spatially 
variable inclusion phase(s). Despite its demonstrated utility, however, the original formulation 
of the higher-order theory is computationally intensive. Herein, an efficient reformulation of the 
original higher-order theory for two-dimensional elastic problems is developed and validated. 
The use of the local-global conductivity and local-global stiffness matrix approach is made in 
order to reduce the number of equations involved. In this approach, surface-averaged quantities 
are the primary variables which replace volume- averaged quantities employed in the original for- 
mulation. The reformulation decreases the size of the global conductivity and stiffness matrices 
by approximately sixty percent. Various thermal, mechanical, and combined thermomechanical 
problems are analyzed in order to validate the accuracy of the reformulated theory through 
comparison with analytical and finite-element solutions. The presented results illustrate the 
efficiency of the reformulation and its advantages in analyzing functionally graded materials. 


1 Introduction 

Functionally graded materials (FGMs) are a new generation of composites wherein the microstruc- 
tural details are spatially varied through nonuniform distribution of the reinforcement phase(s), by 
using reinforcement with different properties, sizes and shapes, as well as by interchanging the roles 
of reinforcement and matrix phases in a continuous manner. The result is a microstructure that 
produces continuously changing thermal and mechanical properties at the macroscopic or contin- 
uum level. This new concept of engineering the material’s microstructure marks a new paradigm in 
both the materials science and mechanics of materials areas since it allows to fully integrate both 
the material and structural considerations into the final design of structural components. 
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Although the area of FGMs is not very old, an enormous amount of research has already taken 
place and the field continues to expand rapidly. Many papers on different aspects of FGMs have 
been published in the past decade in regular journals, special issues of journals devoted to FGMs, 
conference proceedings, and monographs, cf. Pindera et al. (1994a, 1995a, 1997a), Needleman 
and Suresh (1996), Ilschner and Cherradi (1994), Shiota and Miyamoto (1997), and Suresh and 
Mortensen (1998). Therefore, a comprehensive review of the different research activities is outside 
the scope of this report. In keeping with the focus of the undertaken study, an overview is provided 
below of the different approaches employed to model the thermomechanical response of FGMs. 

The analysis of functionally graded materials is a difficult task because of the arbitrary varia- 
tion of material microstructure. In order to analyze the response of functionally graded materials 
under given loading conditions, two distinct approaches have been used to date. The first is the 
uncoupled micro-macrostructural approach. In this approach, FGMs are analyzed directly at the 
continuum level, which allows to reduce a given boundary value problem to a system of differen- 
tial equations with variable coefficients. These variable coefficients represent macroscopic material 
property variations expressed as functions of position that are obtained by local homogenization of 
the microstructure or sometimes taken to have specific functional forms. In some cases, the material 
properties are assumed to be piecewise uniform through appropriate microstructural discretization, 
and then the governing differential equations are solved for each layer subject to interfacial conti- 
nuity and boundary conditions. 

There are various micromechanical models that have been used for homogenizing the microstruc- 
ture of functionally graded materials. These micromechanical models include Voigt and Reuss es- 
timates, Mori-Tanaka method, Composite Cylinder Assemblage (CCA) model, Method of Cells, 
etc., and have been described in detail by Aboudi (1991). The uncoupled analysis of functionally 
graded materials involves the determination of effective (macroscopic) properties at a continuum 
point based on the chosen micromechanics model, a step called local homogenization, which are 
then used in the overall analysis of the structure to determine the macroscopic field quantities. 
This two-step procedure essentially neglects the interaction between non-uniformly distributed in- 
clusions and decouples the locally produced effects of microstructural variation. This can often 
lead to potentially erroneous results, especially when the size of the inclusion phase is large with 
respect to the overall dimensions of the composite, the field gradients are high, or when the num- 
ber of inclusions is small. Hence, these models can only be used in limited circumstances such as 
when the size of the inclusion phase is very small in comparison to the overall size of the analyzed 
structural component and the total number of inclusions is large, Pindera et al. (1995). 

Using the uncoupled approach, various researchers have developed different analytical tech- 
niques for studying particular types of problems in functionally graded materials. In general, these 
techniques include either a specific type of loading, a specific geometry or specific variation of ma- 
terial properties. The type of problems studied involve thermal barrier coating and joint problems, 
crack problems, and design and optimization problems. Representative papers in these areas are 
discussed below. 

One of the most important applications of functionally graded materials is in the field of thermal 
barrier coatings. Internal stresses can be reduced and fracture toughness enhanced with appropriate 
spatial variation of ceramic and metallic phases. Approaches employed to study temperature fields 
and resulting stresses in thermal barrier coatings depend on whether elastic or inelastic analyses are 
conducted, as discussed recently in an extensive review paper by Noda (1999). In the case of inelas- 
tic effects, approaches range from one-dimensional finite-difference analyses, Kokini and Choules 
(1995), to two-dimensional or axisymmetric finite-element analyses based on layer- wise discretiza- 
tion of the coating’s microstructure with piece- wise varying properties, Jian et al. (1995), Delfosse 
et al. (1997). The properties of the individual layers are typically obtained using either simple 
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rule-of-mixtures or Mori-Tanaka estimates. In the absence of inelastic effects, analytical techniques 
and their numerical implementations have been employed with the property variation described by 
continuous functions. For instance, Jeon et al. (1997) discuss an axisymmetric thermoelastic solu- 
tion for an inhomogeneous material by introducing a thermoelastic displacement potential function. 
Sutradhar et al. (2001) use the Green’s function approach in a boundary element setting in order to 
study the 3-D transient heat conduction in functionally graded materials with exponential thermal 
conductivity variation. 

Analysis of graded joints has also received considerable attention due to the large interlaminar 
stresses that may arise along a bi-material interface at the free edge. For example, Drake et al. 
(1993) and Williamson et al. (1993) employed the finite-element method to study the residual 
stresses that develop at graded ceramic-metal interfaces joining cylindrical bodies made of metallic 
and ceramic components. The gradation was modeled using a series of perfectly bonded cylindrical 
layers, with each layer having slightly different properties. The authors’ results demonstrated the 
importance of plasticity effects in the analysis of graded and non-graded interfaces. The authors 
also showed that, in some cases, optimization of the microstructure of graded layers is required to 
achieve reduction in certain critical stress components that control interfacial failure. Along similar 
lines, Suresh et al. (1994) studied the response of elastoplastic bimaterial strips subjected to cyclic 
temperature variations. The authors showed that plastic flow along the material interface at the 
free edge can be modified substantially by altering the constraints at the strip edge. 

Functionally graded microstructures can be very useful in enhancing material’s fracture resis- 
tance through the mechanism of crack blunting and crack-path deflection. In fact, the enhanced 
fracture toughness of graded coatings is due to the local material heterogeneity. In order to keep 
the problem tractable, local material heterogeneity is replaced by spatially varying homogenized 
properties, reducing the problem of a crack in an heterogeneous material to an inhomogeneous 
medium crack problem. Crack problems in inhomogeneous materials have been studied by a num- 
ber of people since the 1960’s. However, systematic studies of this class of crack problems were 
initiated in the early 1980’s by Erdogan and co-workers, cf., Delale and Erdogan (1983), Erdogan 
(1985), Erdogan et al. (1991), and subsequently extended to functionally graded materials. In 
applications involving TBCs, different types of cracks may arise, including cracks perpendicular 
and parallel to the coating’s surface and at the free edge between the substrate and the graded or 
layered coating. These different scenarios have been discussed by Erdogan (1995), Kadioglu and 
Erdogan (1995), Erdogan and Wu (1996), Schulze and Erdogan (1998), Lee and Erdogan (1998), 
and Jin and Paulino (2001). Additional solutions to crack problems in the presence of thermal and 
mechanical loading have been provided by Noda and Jin (1993) who discuss thermal stress inten- 
sity factors for a FGM strip using Fourier transforms, showing the efficiency of suitable material 
selection in reducing these intensity factors. Crack problems in viscoelastic functionally graded 
materials have recently been addressed by Paulino and Jin (2001). 

The ultimate objective in designing a structure for practical applications is not just the ability 
to analyze a given heterogeneous component, but the identification of an optimum design which 
produces the best stress distribution for the given application. Various researchers have provided 
analytical solutions for specific types of problems which are useful in optimizing the design of 
graded structures. These analytical solutions also serve as benchmark solutions for the validation 
and verification of finite-element and other approximate techniques which have been developed or 
are being developed for the analysis of functionally graded materials. For instance, Salzar and 
Barton (1994) incorporated the analytical solution for the axisymmetric elastoplastic response of a 
multilayered cylinder, developed by Pindera et al. (1993), into a commercial optimization code to 
minimize residual stresses in metal-matrix composites using graded interfaces. Horgan and Chan 
(1998, 1999) developed solutions for problems involving pressurized hollow cylinders, rotating disks, 
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and bars under torsion with continuously graded isotropic linearly elastic materials. Ootao et al. 
(1998) discuss the optimization of material gradation for a hollow circular cylinder under thermal 
loading. They use the concept of a laminated composite cylinder with piece-wise homogenous and 
uniform material properties together with the neural network approach in order to determine an 
optimum design. Nadeau and Ferrari (1999) discuss microstructural optimization of a functionally 
graded transversely isotropic layer. Noda (1999) identifies optimal composition profiles to reduce 
thermal stresses in FGMs using perturbation methods. 

The second approach employed in the analysis of functionally graded materials is the coupled 
approach. In this approach, the effects of microstructural variation and the interaction between 
nonuniformly distributed inclusions are explicitly taken into account in the course of solving the 
governing differential equations. This in turn, makes possible the analysis of heterogeneous mate- 
rials with different microstructural scales, in contrast to the uncoupled approach which is limited 
to materials with very fine microstructures. 

The higher-order theory for functionally graded materials (HOTFGM) is an approximate cou- 
pled approach based on a particular volume discretization of the material’s microstructure and an 
averaging approach in the solution of the governing held equations in each subvolume. This theory 
has been developed in a sequence of papers dating back to 1993, Aboudi et al. (1993), in order 
to circumvent the limitations of the uncoupled approach. Summaries of the different stages of the 
higher-order theory’s development have been provided by Pindera et al. (1995b; 1998), and most 
recently in a comprehensive article by Aboudi et al. (1999). As discussed in this review article, 
the theory has been employed to analyze a number of technologically important problems ranging 
from thermally-induced free-edge interlaminar stresses in cross-ply laminates, optimization of fiber 
spacing in laminates subjected to thermal gradients, and thermal barrier coatings. The focus of 
these applications was the demonstration of microstructural coupling effects in functionally graded 
materials as a function of the microstructural length scale, and when these effects can be neglected, 
as well as the demonstration of the utility of functionally graded microstructures in enhancing the 
performance of plate-like structural components subjected to through-thickness thermal gradients. 
Specifically, the higher-order theory has been applied to the following technologically important 
problems: 

• Investigation of the effect of microstructure on thermal and stress fields in MMC plates and 
cylinders 

• Investigation of the use of functionally graded architectures in reducing edge effects in MMC 
plates 

• Optimization of functionally graded microstructures in MMC plates and cylinders 

• Development of guidelines for the design of special coatings in exhaust nozzle applications 
under NASA/Pratt & Whitney Space Act Agreement 

• Investigation of microstructural effects in functionally graded TBCs 

The results obtained so far have demonstrated that this theory is an accurate, cf. Pindera and 
Dunn (1997), and easily implementable tool in the analysis and design of FGMs. Furthermore, 
comparison of the results obtained from the standard micromechanics approach with those of 
HOTFGM has demonstrated the need for a theory like HOTFGM, which explicitly couples the 
micro (local) and macro (global) effects in the analysis, Pindera et al. (1994b; 1995c). 

The recent developments of HOTFGM include extension to cylindrical ccordinates to enable 
analysis, design, and optimization of structural components found in aircraft engine applications, 
Pindera and Aboudi (2000). 
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1.1 Objectives and Outline of the Completed Investigation 

As demonstrated through its applications, the higher-order theory is an accurate and useful ap- 
proach for the analysis of functionally graded materials which takes the micro-macrostructural 
coupling into account and does not neglect the effects of microstructural variation. However, the 
higher-order theory is computationally intensive in those cases which require detailed volume dis- 
cretization in order to mimic realistic microstructural details in certain types of functionally graded 
materials or to capture very high thermal and stress gradients. Such detailed modeling leads to 
a large number of equations whose solution becomes computationally prohibitively expensive and 
sometimes impossible to execute due to the large computer storage requirements. 

This report describes an efficient reformulation of the higher-order theory which leads to a sig- 
nificant reduction in the number of equations required in a problem’s solution. This reformulation 
is accomplished by making use of the so-called local-global conductivity and stiffness matrix ap- 
proaches employed in conjunction with a simplied manner of volume discretization developed by 
Zhong and Pindera ( 2002 ). The simplified volume discretization provides the basis for the refor- 
mulation by enabling the derivation of closed-form relations between surface- average heat fluxes 
and tractions and the corresponding temperatures and displacements associated with a generic 
subvolume. Therefore, volume- averaged quantities employed in the original higher-order theory 
are replaced by surface- averaged quantities as the fundamental unknowns in order to employ the 
local-global conductivity and stiffness matrix approaches in the reformulation. 

The original two-dimensional formulation of HOTFGM is described briefly in Section 2 . The 
motivation for efficient reformulation of the higher-order theory and the reformulation approach are 
then discussed in Section 3 . Section 4 describes the effect of mesh discretization on the field variables 
for selected loading conditions and also provides validation of the reformulated higher-order theory 
by considering several test cases, including the classic Eshelby problem of a circular inclusion 
in an infinite matrix subjected to uniform far-field loading. The results are compared with the 
analytical solution and a finite-element solution obtained using ANSYS. Section 5 demonstrates the 
usefulness of the higher-order theory in the analysis of functionally graded materials by considering 
a thermal barrier coating application. Section 6 summarizes the present accomplishments and 
provides suggestions for future work which should be pursued in this area. 

2 Higher-Order Theory: Original Formulation 

The version of the higher-order theory for materials functionally graded in two directions, or 
HOTFGM- 2 D, is based on a geometric model of a heterogeneous composite graded in the X2 — %3 
plane which occupies the region 0 < X2 < H, and 0 < x% < L, Fig. 1 . The microstructural 
pattern is represented by discretizing the cross-section of the heterogeneous composite into N q and 
N r generic cells in the intervals 0 < X2 < H, and 0 < x% < L, respectively. The indices q and r 
(q = 1 , 2, ..., N q and r = 1,2, ..., 7V r ), identify the generic cells in the X2 — X3 plane. The generic cell 
(<7, r) consists of four subcells designated by (/Jy), where each index /?, 7 takes a value 1 or 2 which 
indicates the relative position of the given subcell within the generic cell along the X2 and X3 axis, 
respectively. The thermomechanical properties of the material within a subcell are assumed to be 
constant. The composite is assumed to be infinite in the x\ direction, whereas the dimensions of 

(q) (r) (V) 

the generic cell along the functionally graded directions X2 and x% are /q , , and •> and 

can vary arbitrarily such that 

N q N r 

= + L = (1) 

q= 1 r= 1 
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Fig. 1. H0TFGM-2D representation of a composite functionally graded in the X 2 and X 3 
directions with uniform microstructure in the x\ direction. 


The composite is subjected to combined thermomechanical loading on its bounding surfaces in 
the X 2 — X 3 plane given in terms of temperature or heat flux, and surface displacements and/or 
tractions. Further, the strain e\\ (Eu = 0, for plane strain) is uniform in the x\ direction. Given the 
applied thermomechanical loading, an approximate solution for the temperature and displacement 
fields is constructed. The solution strategy involves volumetric averaging of the field equations 
together with the imposition of boundary and continuity conditions in an average sense between 
the sub-volumes used to characterize the material’s microstructure. As described in Subsections 
2.1 and 2 . 2 , the temperature and in-plane displacement fields in each subcell of a generic cell 
are approximated using a quadratic expansion in the local coordinates x^ \ placed at the 
subcell’s centroid. This temperature and displacement field representation is sufficient to capture 
the local effects created by the thermomechanical field gradients and the microstructure of the 
graded material with finite dimensions in the functionally graded directions. 


2.1 Thermal Analysis 

The thermal boundary conditions specified in terms of temperature or heat flux distribution on the 
bounding surface in the X 2 —X 3 plane are assumed to be steady-state. The temperature distribution 
in the subcell (/?y) of the (< 7 , r) th generic cell measured with respect to a reference temperature T re j, 
is denoted by This temperature field is approximated by a second-order expansion in the 

local coordinates x^\ as follows: 


1 -{P)rp((3 7 ) , -('y)rri(Pj) , ^ P \rn(Pl) , ^ /o -G ) 2 _ ^7 1 

1 ~ 1 (m) X 2 1 (10) ^ X 3 1 (01) + 0*^2 A ) 1 (20) A ) 1 ( 




i(r) 2 


(0 7 ) 

-( 00 ) ^^2 -*■ ( 10 ) “^3 ( 01 ) ^ 2 V °^ 2 4 ; J _ ( 20 ) 4 7 ^ ( 02 ) 


(2) 
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where is the volume-averaged temperature in the subcell (/?y), and T (mn) (™,n = 0,1, or 

2 with m + n ^ 2) are the higher-order coefficients which provide a better approximation of the 
temperature held throughout the subcell. 

For a generic cell with four subcells containing arbitrarily specified materials, five unknowns 
(i.e., are associated with each subcell, producing twenty unknowns for each generic 

cell. This results in 20N q N r unknown quantities for a composite with N r rows and N q columns of 
generic cells. These unknowns are determined by first satisfying the zeroth, first and second moment 
of the heat conduction equation in each subcell in volumetric sense. Subsequently, continuity of 
heat flux and temperature is imposed in an average sense at the interfaces separating adjacent 
subcells, as well as neighboring cells, followed by the boundary conditions. 


2.1.1 Heat conduction equation 


Under steady-state thermal boundary conditions in the X 2 — £3 plane, the heat flux held in the 


material occupying the subcell (/?y) of the (g, r) th generic cell, in the region 
7 , must satisfy the steady-state heat conduction equation given by 


Jfi) 


x , 


-2 ri (3 > 


X 


( 7 ) 


< 


B4 0l) , 

dx 2^ dx ^ 


= 0 , (/?, 7 = 1 , 2 ) 


(3) 


The heat hux components at any point passing through a subcell (/?y) are derived from the 
temperature held according to the Fourier’s law of heat conduction given by 


= _ k m 


QT (£7) 

dxf 1 


(i = 2, 3; no sum ) 


(4) 


where are the heat conductivity coefficients of the material in the subcell (/?y) assumed to be 
orthotropic, and no summation is implied by repeated Greek letters in the above and henceforth. 

In order to satisfy the above steady-state heat conduction equation in a volumetric sense, the 
following volume- averaged quantities are dehned in the higher-order theory 


nO ^7) = 


h^/2 ,jyv 2 


7) , 


4 (<?>'•) 


-hf /2 J •-!%> /2 


(r) 


[x^) m (x^ ) n q^^ dx 9^ dx i 7 ^ 


(5) 


where is the area of the subcell (/5y) in the (g, r) th generic cell. ^ is the average 

value of the heat hux component qf ^ in the subcell (/?y), whereas for the other values of m and 
n Eq. (5) dehnes higher-order heat huxes. Satisfaction of the heat conduction equation given by 
Eq. (3) results in the following four conditions in terms of the volume- averaged hux quantities 


Q ( 2%/^ + Q'm i>A 


( 07 ) 


(q, r ) 


= 0 


(6) 


Equation (6) can be expressed in terms of the coefficients by evaluating the hux quantities 

dehned in Eq. (5) . Substituting Eq. (4) into Eq. (5) , and performing the required volume 
integration yields the following non-vanishing zeroth-order and hrst-order heat huxes in terms of 
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the coefficients T^\ 

( mn ) 


Q 


Q 


i(07) 

u(M)rp(Pi) 

* 1 (10) 

h(l) 2 

lSPl) (3 rpiP'l) 
K 2 4 J (20) 

2(0,0) 

i(07) _ 

2(1,0) 

i(07 ) 

^3 J (oi) 

,(r)2 

U (Pi) 7 rp(fry) 
3 4 (02) 

3(0,0) 

\(Pi) _ 

3(0,1) 


(7) 

( 8 ) 
(9) 

( 10 ) 


2.1.2 Heat flux continuity conditions 

The heat flux continuity at the interfaces separating adjacent subcells within the generic cell (g, r) 
is imposed in an integral sense as follows 


1 


r4 r) / 2 




4 17) (^.4 7| 1 


(lx) 


dx o 7 ^ = -fy 

6 ,(r) 


1 


rW/ 2 

'-W/2 


% 


(27)/ 


h 


(<?) 




->®3 


(?,r) 


cfcc 


~(7) 


A q) /2 

hW J-hf/2 


1 

y 


(r) l(9,r) 

,f )(4«, !l_) I 


= 


/i (g) 

n P 


A q) /2 

l-h^/2 


% 


092) , 09) J 

3 l x 2 > 


(r) 


(9> r ) 


dx, 


( 13 ) 


( 11 ) 

(12) 


Similarly, the heat flux continuity at the interfaces between generic cells is satisfied in an integral 
sense by requiring that 


1 

rW/ 2 

/M 

/-4 r) /2 

i 

r^ 9) /2 

/,(«) , 
"73 

/-4 9) /2 


4 l7) (-^ 


hl ? +1) Jyh 

5 x 3 v 


-| (9+1, r) 


dx i 7 ^ = 


r4 r) /2 


^3 V x 2 ’ 


(r+1) 1 


= 


y-jW/2 

r h f/ 2 

hf J-hf / 2 


% 


( 27 ) 


h. 


(9) 


(7)> 


1 («)»•) 


■^3 


dx^ 


m (f (P) k_ 

% \ x 2 1 o 


(r) 1 («> r ) 


dx, 


=08) 


(13) 

(14) 


Satisfaction of the above heat flux continuity conditions (11) — (14) results in the following eight 
equations in terms of the volume-averaged flux quantities 


12^2(1, 0)/^ 11 + ^2(0,0) ^2(1, 0)/^ 2 

(9,r) 

_ /!,„ | 

2(0,0) ^ 2^ 2 (°>°) 


427 ) 


427 ) 


1 (ix) 


^2(0,0) ^(iV^ 2 




= 0 


1 

+ 2 


^2(0,0) + 6( 52(1,0)/^ 2 
n(^2) , «n (/32 ) //„ 

h?3(o,0) Dt » 3(0,1) / t2 

q(/32) . n/-^(l32) j 7 

“*3(0,0) D ^3(0,l)/ t2 


l(27) 


12^3(0,1)/^ 9“ ^3(0,0) ^3(0,1)/^ 


l082) 


(9,r) 


1 d_/q(/^2) 

h! 3(0,0) -r 2 ^3(0,0) 


^3(0,1) /^ 2 


1 (?, r ) ^ 

+ X 


(9-l.r) 


= 0 


(g.r-1) 


= 0 


(9>r-l) 


= 0 


(15) 

(16) 

(17) 

(18) 


Equations (6) and (15) — (18) have been obtained in terms of the volume-averaged flux quantities 
after some complex algebraic manipulations. For detailed derivation of these equations, refer to 
Aboudi et al. (1996) . These equations can now be expressed explicitly in terms of the microvariables 
^(mn) the use of Eqs. (7) — (10). 
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2.1.3 Temperature continuity conditions 


The temperature continuity at the interfaces separating adjacent subcells within the generic cell 
(q, r ) is imposed in an integral sense as follows 


f4 r) / 2 


4 r) 


/2 


7,(9) . . 

r^)(V,4 7) ) 


(ix) 


dX = —7~T 


r-4 r) /2 


4 r) y-d r) 


/2 


r( 27)(_!!2_ 


/,(«) / X 

) x 3 J 


(9,r) 


da: 


(7) 




^ Th^/2 




(r) 1 (?> r ) 


da; 


09) _ 


^ } /2 


2 ’(«). /_».(«) , 


h 


P 


'-hr/2 


r (/3 2 )( S (/y_k-) 


(r) 1 (?> r ) 


da: 


~(0) 


(19) 


( 20 ) 


Similarly, the temperature continuity at the interfaces between generic cells is satisfied in an integral 
sense by requiring that 


1 rW/ 2 


J'( 1 7)^_ 


/l 


(9+1) 


1 


r^ 9) /2 


^ 9) -z-^/ 2 


2 

x 2 


™(t) > 

• x 3 ' 


(9+1, 0 1 „/W / 2 

= J_ / 7 7 
3 J-ip/2 


/,(<?) 


(qx) 


dx . 


~( 7 ) 


y(/9i)^(^) ‘l 


(r+l) 1 (9,r+l) 


da: 


~(0) - 


-h^/2 


2 ' 0) ./_*.(«) , 


h 


P 


'-hr/ 2 




if) .l (, ' r) 


dx, 


03) 


( 21 ) 


(22) 


Satisfaction of the above thermal continuity conditions results in the following eight equations given 


explicitly in terms of the microvariables T, 


(07) 


(ran) 


7^(17) 1 U 7-1(17) 1 1 7,27-1(17) 
J (oo) ^ 2 AtlJ (i°) ^ 4 n i J ( 20 ) 


7,(27) . 7-1(27) . ^,27-1(27) 

J (00) "l" 2 ^ J (10) 4 ^ 1(20) 


7i(/31) 1 jj_/ rp(Pl) 1 I,2 T (/31) 
J (oo) ^ 2 tlJ (oi) _l " 4 ,, 1 J ( 02 ) 

r r(P2) 1 1, r 092) 1 \i2rp{p2) 
1 (00) ^ 2 t2J (01) _l " 4 t2 -l(02) 


- («,»•) 
(9,r) 
(9,r) 
(9.’’) 


p(27) 

( 00 ) 


^ /, „7 1 ( 2 7) 1 ^ ^ 27 ,( 27 ) 

2^ 2J (10) + 4 n 2 l( 20 ) 


1 (9,r) 


p(l7) _ U 7 - 1 ( 17 ) 1 17 , 27 ,( 17 ) 
-(oo) 2 ,il_t (i 0 ) 4 At i x (20) 


t; 


(d2) 


( 00 ) 


7-1 (/31) 
( 00 ) 


1 7 , rp(p2) 1 -( l2rp{P2) 

' 2^ 2i (01) + 4 <2 ( 02 ) 

t/,7-'(dl) 1 1/27,031) 
2 t i J (oi) 4 t l-l(02) 


- (g+l,r) 

- (9,r) 

(9, r+l) 


(23) 

(24) 

(25) 

(26) 


2.1.4 Boundary conditions 

Some of the heat flux and temperature continuity conditions described above are not valid for cells 
located on the boundaries defined by the indices q = 1, A^, and r = 1, 7V r . For the set of boundary 
cells with q = 1, heat flux continuity between the given generic cell and preceding generic cell, 
Eqs. (15) — (16), is not applicable. Similarly for q = N q , temperature continuity between the given 
generic cell and following generic cell, Eq. (24), is not applicable. These conditions are replaced 
by the continuity of heat flux between adjacent subcells within the generic cell (l,r), and in case 
of thermal boundary conditions by 
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where T left( x 3 )> and T^ghJxz) are the piece- wise uniform applied temperatures on the external 
boundaries in the X 2 — X 3 plane. Similar reasoning holds for generic cells (q, 1), and ( q,N r ). If 
the temperatures are defined at the top and bottom surfaces for instance, the applied boundary 
conditions are given by 


h 


(?) 


h 


(?) 


A q) / 2 
'-hf/2 

r ^/2 

l-W /2 


J'(h 1 )( 


f (P) 


/(l) 


(?.l) 


dx 


JP) 


= T, 


(?) 


bottom 


( x 2 ) 




(N r ) 1 i^ N r) 


dx 


™(0) r r(. < i) 


= 


(29) 


(30) 


If the heat flux at any of the boundaries is defined instead of the applied temperature, Eqs. (27) — 
(30) are modified accordingly. 


2 . 1.5 Solution for T^l 

(mn) 

To obtain a solution for the 20N q N r unknown coefficients (T^^ in each (/?y) subcell), 20N q N r 
equations are required. These equations are assembled using the governing field equations, Eq. 
(6), heat flux continuity conditions, Eqs. (15) — (18) , and the temperature continuity conditions, 
Eqs. (23) — (26) , together with the imposed boundary conditions on the external boundaries of the 
composite, Eqs. (27) — (30). The final system of equations obtained is symbolically written as 

k T = t (31) 


In the above equation, k is the structural thermal conductivity matrix which contains information 
on the geometry and thermal conductivities of the individual subcells (/?y) in the N q N r cells. The 
thermal coefficient vector T contains the unknown coefficients that describe the temperature field 
in each subcell, i.e., 


T = 


.(ii) 

■li 


T (22) i 

*•*’ ± N q N r j 


where 

Tg 7) = [r ( 00),T (1 0),r (01) ,T (20) ,T ( 0 2) ]^ 


The thermal force vector t contains information on the boundary conditions. 


2.2 Mechanical Analysis 

The next step is to find the displacement and stress fields in the heterogeneous composite due to the 
temperature field obtained in Subsection 2.1, and/or mechanical loading at the external boundaries 
applied in a manner which is consistent with the global equilibrium requirements. Towards this 
end, the displacement field in the subcell (/3y) of the (#, r) th generic cell is approximated by a 
second-order expansion in the local coordinates xip and as follows: 

uf l] = x\E\\ (32) 


u , 


(07) 


(07) j_^(0)u/(07) _l ^(7)pjy/07) 


kk 2(oo) ^ x : 


^2(10) + X 3 


2 ( 01 ) 


+ - 


,(9)2 


-)W® + i(34 7,z - 


i ;( r ) 2 

1 h \uAPi) 


)^2(02) (33) 


(9)2 


U (P A — + x^W^^ + ~(3x^ 2 — — -I- ~(3x^ 2 — — )W 

u 3 - *^3(00) + x 2 VV 3(W) + X 3 VV 3(01) + o ^*2 „ ) VV 3(20) + o { 6X 3 a > W l 


3(00) 


3(10) 


3(01) 


3(20) 


{Pi) 

'3(02) 


(34) 


NASA/CR— 2002-21 1909 


10 



are 


where W^ggj are the volume-averaged displacements in the subcell (/3y), and (i = 2,3) 

the higher-order quantities which determine the displacement and stress field at specific locations 
within the subcell. Note that, in Eq. (32) , x\ is associated with the global coordinate system 
fixed to the edges of the composite plate and not the local coordinate system associated with each 
subcell. Also note that there are no product terms of the form x^x^ appearing in the above 
displacement field approximation due to the employed volumetric and surface- averaging procedure 
technique. The quantity £n is the uniform strain in the x\ direction which is zero for plane strain 
analysis. For generalized plane strain, e\i is determined from the condition 



a, 


(07) 






= 0 


(35) 


The other 10 unknown coefficients W^n) associated with each subcell (/?y) of the (< 7 , r) th generic 
cell are determined from the conditions analogous to those employed in the thermal problem. Here, 
the heat conduction equation is replaced by the two equilibrium equations, and the conditions 
involving continuity of heat fluxes and temperature at the interfaces are replaced by the continuity 
of tractions and displacements. 


2.2.1 Equilibrium equations 

The stress field in the subcell (/?y) of the (< 7 , r) th generic cell generated by the given temperature 
field and the applied mechanical loading must satisfy the equilibrium equations 


Q a (07) q (07) 

acje 2 j aa 3j 

dx 2 ^ dx ^ 


0, O' = 2, 3) 


(36) 


For an orthotropic elastic material occupying the subcell (/3y) of the (< 7 , r) th generic cell, the stress 
components are related to the strain components through the familiar Hooke’s law 


(07) _ >o(07) J07) _ T(J3y) 
a ij ~ ° ijkl S kl a ij 


(37) 


where Cijki are the stiffness tensor elements of the material in the subcell (/Ty), are the elastic 
strain components, and are the components of the so-called thermal stress given by 


T(07) _ p(07)^(0v) 
ij ij 


(38) 


where rO 7 ^ is the product of the stiffness tensor components and thermal expansion coefficients. 
7) represents the temperature change at a particular location in the subcell. 

The components of the strain tensor in the individual subcells are obtained from the following 
strain-displacement relations 


.(07) _ 


1 du 


(07) 


-ij 


dx, 


*(■) 


+ 


du i 

dx 


(07) 


0 ) 


(39) 


As in the thermal formulation case, in order to satisfy the equilibrium equations in an average sense 
the following volume- averaged stress quantities are defined in the higher-order theory 


g (07) 

0 ij(m,n ) 


hf/ 2 /< 72 


(r) 


A 


(q, r ) 

(^ 7 ) 


-h { g q) /2J-l\ r, /2 


( r ) 


( 4 0) ) 


(x^Ta^dx^dx^ 


(40) 
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where 0) is the average value of the stress component in the subcell, whereas for other 

values of m and n Eq. (40) defines higher-order stresses that are needed to describe the governing 
field equations of the higher-order continuum. Satisfaction of the equilibrium equation (36) results 
in the following eight equations in terms of the volume-averaged stress quantities 


q(Pi) 

3 2i(l,0 


)/h 2 p + 


s: 


(Pi) 

3.? (0,1) 


A 


(ix) 


= 0 


(41) 


Equation (41) can be expressed in terms of the coefficients by evaluating the stress quan- 

tities defined in Eq. (40) . Substituting Eqs. (37) , (38) , and (39) in Eq. (40) , making use of the 
displacement field representation, Eqs. (32) — (34), and performing the required volume integra- 
tion yields the following non-vanishing zeroth and first-order stress components in terms of the 
coefficients 


c(/E) 

°22(0,0) 

- °22 |/|/ 2(10) ^ °23 ^3(01) 1 22 1 (00) 

(42) 

qOE) 

°22(1,0) 

1 rp(Pl) IJt(Pi) ^ L(^)2-p(/37)ryn(/37) 

- 4 a /3 °22 ^2(20) i2 n P 22 ( 10 ) 

(43) 

o(/E) 

°22(0,1) 

^ ^ j(r)2'p(P'y)rp(P'y) 

- 4 S °23 ^*2(02) 12 7 22 ( 01 ) 

(44) 

q(/E) 

°33(0,0) 

- °23 ^2(10) +°33 ^3(01) 1 33 1 (00) 

(45) 

qOE) 

°33(1,0) 

_ ^ AqA rp{dl) Tjr(dl) ^ ij(9)2-p(^7)/yi(^7) 

- 4 ^/3 °23 ^2(20) 12 n & 33 1 (10) 

(46) 

qOE) 

°33(0,1) 

_ ^ ]{r)2 1 /(r)2-p(^7)a-i(^7) 

- 4 S °33 ^2(02) 12 1 33 ( 01 ) 

(47) 

e(/E) 

°23(0,0) 

_ _L w(^7) ^ 

- 0 44 1^2(01) + ^(lO)^ 

(48) 

q(P 7) 
°23(1,0) 

__ 1 7 (9)2 r y(^7) T J/(/^7) 

- 4 a /3 °44 ^3(20) 

(49) 

qOE) 

°23(0,1) 

__ 1 i{r)2 r (Pl) w (Pl) 

- 4 S °44 ^2(02) 

(50) 


2.2.2 Traction continuity conditions 

The traction continuity at the interfaces separating adjacent subcells within the generic cell (g, r) 
is imposed in an integral sense as follows 


1 r4 r) /2 

W ^-4 r) /2 

A 9) / 2 

hf J-hf ! 2 


7 (q) 

_(E )( [h Ai) 

a 2j v 2 ’ x 3 


l («>»•) 


(/31) f (/3) 
a i3 \ x 2 > o 


(W,l ( " r) 


dx^ 

dx^ 


1 r4 r) /2 

W 

hf J-hf/2 


(27)/_^2_ ™ 
°2 j V 2 ’ X 3 


k. 


(?) 




-| (q,r) 


dx 




a. 


m ( =09) ^2 

V x 2 ’ o 


(r) l(^ r ) 


z3 


dx . 


09) 


(51) 

(52) 


where i, j = 2,3. Similarly, the traction continuity at the interfaces between generic cells is satisfied 
in an integral sense by requiring that 


/O) 


X V) / 2 

I -1^/2 


a 


7 (<7+l) 

(EE "i ^(t) > 


-1 (g+l,r) 


21 




dXo 7 ^ = -pr 
3 7 ( r ) 


^E/2 

-4 r) /2 


CT 


Aq) 

(27) / ? '2 
21 1 9 ' 




1 Ed) 


da: 


( 7 ) 


(53) 
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r h P q) / 2 

hf J-hf/2 


a, 


(01) /=(/*) 


i3 


(®2 


7 


(r+ 1 ) 


(«,»•+ 1) 


dx 


Ad) 


j_ r h f/ 2 
hf J-hff 2 


cr, 


(02), (0) l A 

i3 ( x 2 s r 


(r) 


(qx) 


dx 


Ad) 


(54) 


where i,j = 2, 3. Satisfaction of the above traction continuity conditions (51) — (54) results in the 
following sixteen equations in terms of the volume-averaged stress quantities 


O )/^ 1 + ^ 2 j ( 0 , 0 ) ^^ 2 i ( l , 0)/^ 2 

n (g>»0 


,(27) 


427) 


(9,4 


S^^ + 65^^7/12 


427) 

J 2j( 0,0) 


4j(l,0)/ 


(9-1, r) 


= 0 


_ c (17) 1 c( 2 7) _ qcW /u 

°2j(0,0) + 2 2j(0,0) ^(l.O)/" 2 

_-|oq(01) n 1 <40 2) _ «<40 2 ) / 7 

iZ,5 3j(0,l)/ (l _l " D 3j(0,0) D '- , 3j(0,l)/ t2 


1 

+ 2 


sAL + 6S!%} n ,/h2 


1 (9-l.r) 


'2j(0,0) 


(9,r) 


+ 65,^2 ^/i 2 


4/32) 
^3.7 (0,0) 


1 ^ o(/^2) _oq(/32) / 7 

3J (0,0) ^ 2 3 J'(0,0) OO 3j(0,l)/^ 2 


-sf‘> 


(9,r) 


1 

+ 2 


q(d^) 1 cc(42) M 

D 3j(0,0) ^ 0,3 3iYn.1V J 2 


"2j(l,0)/ 

402) 

5 3i(0,l)/ 

4/32) 

J 3i(0,l)/ 


= 0 


(g,r-l) 


= 0 


1 (9,r-l) 


= 0 


(55) 

(56) 

(57) 

(58) 


where j = 2, 3. Equations (41) and (55) — (58) have been obtained in terms of the volume- averaged 
stress quantities after some complex algebraic manipulations. For detailed derivation of these 
equations, refer to Aboudi et al. (1996) . These equations can now be expressed explicitly in terms 
of the microvariables an d ^s(mn) w ^h the use of Eqs. (42) — (50) . 

2.2.3 Displacement continuity conditions 

The displacement continuity at the interfaces separating adjacent subcells within the generic cell 
(q, r ) is imposed in an integral sense as follows 


r4 r) / 2 


J-W 


/2 


(! 7 ) ( 7 ) 

2 V 2 5 X 3 ) 


u . 


(<?x) 


dx i 7 ^ = — 7— r 


rl ^ /2 


4 r) 7-^ 


/2 


.427)/^ =(7)> 
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^(7) 


A q) A 

J-hf / 2 


(r) 1 («.»•) 
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n p 


'~h ( g q) / 2 


u 


m<Ad) 


ik 


(r) l(9- r ) 


^2 5 
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(59) 

(60) 


Similarly, the displacement continuity at the interfaces between generic cells is satisfied in an 
integral sense by requiring that 


1 

47/2 

/M 

l/'y 

/-i7/2 

1 
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^(«) , 
n f3 

l-h^/2 


h 

.(17) / 'h =(7)> 






4+i) 
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47/2 
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^(7) 




/(r) 

(I®)/ 09) 7\ 

3 V x 2 > 2 ' 


(9,r) 


dx, 


~( 0 ) 


(61) 


(62) 


Satisfaction of the above displacement continuity conditions results in the following sixteen equa- 


tions expressed directly in terms of the microvariables W< 


(07) 


2 (ran) 


and 117 


(07) 


3 (ran) 


114 


(17) 


,(00) + + j^po) 


41,7) 


(9,r) 


tj4 2 7) _ tj/( 27) , I/,2iy4 2 7) 

7(00) 2 A(10) + 4 ^2 7/(20) 
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(63) 
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(64) 
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(65) 

( 66 ) 


2.2.4 Boundary conditions 

As in the thermal analysis case, some of the displacement and traction continuity conditions 
described above are not valid for generic cells located at the boundaries defined by the indices 
q = 1, N q , and r = 1, N r . For the set of boundary cells with q = 1, the traction continuity between 
the given generic cell and preceding generic cell, Eqs. (55) — (56), is not applicable. Similarly for 
q = N q , the displacement continuity between the given generic cell and following generic cell, Eq. 
(64), is not applicable. These conditions are replaced the continuity of tractions between adjacent 
subcells within the generic cell (1, r), and in case of traction boundary conditions by 


4 r) / 2 


4 r) J-W 


/ 2 


4 7) (- 


4 ] Jn) 
2 ’ 3 


1 (!.»•) 


dx i 7 ^ = 


/O) 

i'y 


rWf 2 

I - 1^/2 


' „ ^ h {Nq) 

_(1 l)( n 2 

a 2 j V n 


>4 7) ) 


(N g ,r) 


t{ lefA X * ) 


dx i 7 ^ = t 


O') 

right 


(*3) 


3 = 2,3 


(67) 

( 68 ) 


where t^ t {x^), and t^ ht (x 3) are the piece- wise uniform applied surface tractions. Similar reasoning 
holds for generic cells (g, 1), and (#, N r ). If the top and the bottom surfaces are fixed for instance, 
the applied boundary conditions are given by 


rhf/2 


J-hf ! 2 


rhi q) / 2 


1 

y 

l /3 


u, 


091)^09) Z 1 


(1) “I («•!) 


dx. 


~09) - 




..091)^08) ^2_ 
V x 2 ’ 2 


(JV r ) 1 (^Ar) 


= 0 


(69) 

(70) 


For other type of boundary conditions, Eqs. (67) — (70) are modified accordingly. 


2.2.5 Solution for W^\ 

Thus for the solution of the 40 N q N r unknown coefficients in each (/?y) subcell, 40N q N r 

equations are assembled using the governing field equations, Eq. (41), traction continuity condi- 
tions, Eqs. (55) — (58) , and the displacement continuity conditions, Eqs. (63) — (66) , together with 
the imposed boundary conditions on the external boundaries of the composite, Eqs. (67) — (70). 
The final system of equations obtained is symbolically written as 

K U = f (71) 
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In the above equation, K is the structural stiffness matrix which contains information on the 
geometry and thermomechanical properties of the individual subcells (/3^y) in the N q N r generic 
cells. The displacement coefficient vector U contains the unknown coefficients that describe the 
displacement field in each subcell, i.e., 


u=[<) 


...,u 


( 22 ) 

NqN r 


(72) 


where 

Ug 7) = [W 2m , ...,W 3( 02)]g 7) 

The mechanical force vector f contains information on the boundary conditions and the thermal 
loading effects generated by the applied temperature. 


3 Higher-Order Theory: Efficient Reformulation 

The basic equations of the original higher-order theory outlined in Section 2 will be contrasted and 
compared with the reformulated equations developed in this section. The reformulation simplifies 
the theory and makes it computationally more efficient, as will be demonstrated later. It must 
be emphasized that the basic structure and concepts of the higher-order theory with regard to 
the temperature and displacement field approximation based on the quadratic expansion in the 
local coordinate system attached to a subcell’s center, satisfaction of the governing field equations 
(steady state heat conduction equation and equilibrium equations) in a volumetric sense, and the 
interfacial continuity and boundary conditions in a surface-integral sense, do not change. The 
major changes involve 

• simplification of the volume discretization by eliminating the concept of generic cells, leaving 
only subcells as the basic building blocks of the material’s microstructure 

• replacement of volume-averaged heat flux and stress quantities defined in Section 2 by Eqs. (5) 
and (40) , respectively, by surface-averaged quantities (temperature and heat flux for thermal 
analysis, displacements and stresses for mechanical) associated with subcell interfaces as the 
fundamental unknown quantities 

The above changes set the stage for reformulating the higher-order theory based on the local- 
global conductivity and stiffness matrix approaches described in the following subsections. As 
shown in Fig. 2, the microstructural pattern of the heterogeneous composite functionally graded 
in the X 2 — £3 plane is represented by discretizing the cross-section into Np and AT 7 subcells in the 
intervals 0 < X 2 < if, and 0 < X 3 < L, respectively. The thermomechanical properties within a 
subcell are assumed to be constant. The composite is assumed to be infinite in the x\ direction, 
whereas the dimensions of the subcell along the functionally graded directions X 2 and £3 are hp, 
Z 7 , and can vary arbitrarily such that 


Np tv 7 

H = Y J hp L = (73) 

( 3=1 7=1 

The composite is subjected to a combination of surface displacements and/or tractions applied 
in the X 2 — X 3 plane along with a uniform strain £n ( sn = 0 , for plane strain) in the x\ direction. 
An arbitrary surface temperature or heat flux distribution may also be prescribed. 
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Fig. 2. Reformulated HOTFGM-2D representation of a composite functionally graded in the x 2 
and X 3 directions with uniform microstructure in the x\ direction. 


In the case of the original higher-order theory, the global conductivity matrix relates the heat 
flux defined at the boundaries to the microvariables T^’ 7 ? and the global stiffness matrix relates the 

ymn) ° 

tractions applied at the boundaries to the microvariables W^ ,nf \. In contrast, in the reformulated 
version the global conductivity matrix relates the heat flux defined at the boundaries to the common 
interfacial surface-averaged temperatures and surface- averaged temperatures at the outer surfaces. 
The global stiffness matrix relates the tractions at the boundaries to the common interfacial surface- 
averaged displacements and surface- averaged displacements at the boundaries. 

3.1 Motivation for Reformulation 

The motivation behind the reformulation of the higher-order theory is to decrease the number of 
equations in order to make it computationally more efficient. This is achieved by using the local- 
global conductivity and local-global stiffness matrix approach. This approach involves formulation 
of a local conductivity matrix that relates heat fluxes to temperatures, and a local stiffness matrix 
that relates tractions to displacements, evaluated at the external boundaries of each subvolume, 
and has been described in detail by Pindera (1991) in the context of mechanical boundary-value 
problems. In the reformulated version of the higher-order theory, the local conductivity matrix 
relates the surf ace- averaged fluxes to the surface- averaged temperatures for a particular subcell 
(/?, 7), accounting for the satisfaction of the steady-state heat conduction equation. Similarly, the 
local stiffness matrix relates the surf ace- averaged tractions to the surface- averaged displacements 
for a particular subcell (/?, 7), accounting for the satisfaction of the global equilibrium require- 
ments in a volumetric sense. Once the local conductivity and stiffness matrices are formed, we 
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Fig. 3. Reduction in the size of the global stiffness matrix due to reformulation. 



Fig. 4. Reduction in the size of the global stiffness matrix due to reformulation for equal number 

of subcells in the X 2 and xs directions. 
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use them to construct the global conductivity and stiffness matrices. Here, we enforce the traction 
and displacement (heat flux and temperature for thermal analysis) continuity conditions at the 
interfaces of the adjacent subcells in an average sense, thereby reducing the size of the global 
conductivity matrix and global stiffness matrix by more than fifty percent. 

When the composite is discretized into Np and N 7 subcells in the X 2 and x% directions, respec- 
tively, the rank of the global stiffness matrix in the original higher-order theory is lONpN^ since, 
according to Section 2, each generic cell contains 4 subcells. As we shall see in this section, for the 
same number of subcells in the x 2 and x% directions, the rank of the global stiffness matrix in the 
reformulated version is reduced to 4 NpNry + 2 Np + 2N 1 (ANpN 1 + 2 Np + 2 AT 7 + 1 for generalized 
plane strain). Thus for higher values of Np and AT 7 , the rank of the global stiffness matrix is reduced 
by approximately sixty percent as shown in Figs. 3 and 4. 


3.2 Thermal Reformulation 


The temperature distribution T^’ 7 ^ in the subcell (/?, 7) measured with respect to the reference 
temperature T re f is approximated by the same second order polynomial expansion in the local 
coordinates x^\ x^ as that given by Eq. (2), reproduced below for convenience 

r p{Pil) 7-4/47) 1 ~(P)rp(Pa) I ~(l)rp(Pa) | J^/oA/^) 2 _ h/ 3 \ /o-(7) 2 _ ^7 \rp(Pg) ( 7 A\ 

1 _ ( 00 ) ^ x 2 1 (10) ^( 01 ) ^ 2 4 > 1 {20) 3 4 ' (° 2 ) ^ 

As given in Eq. (4), the heat flux at any point passing through a subcell (/?, 7) is dictated by the 
Fourier’s law of heat conduction, 


9, 


(Pa) 




dx^ 


, (i = 2, 3; no sum) 


(75) 


where are the heat conductivity coefficients of the material in the subcell (/?, 7). Substituting 

Eq. (74) for the assumed temperature field in the above equation and simplifying, we get 


JPa) _ _u(P> 7) ( T (Pa) 1 o-^(P) T (Pa)\ 
Q 2 — K 2 y 1 (10) + ^2 1 (20) ) 

n (Pa) _ 7 ) ( T (Pa) , q T (7) T (/k7)4 

% ~ ^3 W(01) + ^3 (02) ) 


(76) 

(77) 


As mentioned earlier, the reformulation employs surface- averaged quantities in contrast with 
volume- averaged quantities in the original formulation. Hence, the surface-averaged temperatures 
and surface-averaged heat fluxes at the outer faces of the subcell are defined below. The surface- 
averaged heat flux Q 2 going into the left face of the subcell (/?, 7) is defined as 


Q 


-(Pa) 


_ 1 f h/2 „(Pa)( hp (v)x , ( 7 ) 


h J-i^/2 


% 




(78) 


The surface-averaged heat flux going out of the right face of the subcell (/?, 7) is defined as 


o 2 +(/3,7) = 


1 r l i/ 2 
^7 J — l-y/2 


% 


(ffnO ( hp ( 7 )s , ( 7 ) 

, X 3 


(79) 


The surface-averaged heat flux Q 3 going into the bottom face of the subcell (8. 7) is defined as 


«<*” = T P e/2 « t-’hif.-piif' 

h f3 J-hp/2 2 


(80) 
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The surface-averaged heat flux going out of the top face of the subcell (/?, 7) is defined as 


q 3 +(/3,7) 


1 

hp 



l -^)dx { 2 P) 


(81) 


Substituting Eq. (76) into Eqs. (78) and (79), and Eq. (77) into Eqs. (80) and (81) , and performing 
the required integrations, we obtain 


0 

Q 

Q 

Q 


-(A 7 ) _ 

2 “ 

ufin) 

+(^.7) _ 

2 ““ 

ufin) 

-(A7) _ 

3 — 

ju(A7) 

/l, 3 

+(£.7) _ 

3 ““ 

ufin) 

k 3 



_ ^fi rpfi, 7)^ 

2 ( 20 ) ) 



i rpfin)\ 
+ 2 (20) ) 



^7 rp(Pn)\ 

2 (02) J 



, ^jn_rpfi,-y)\ 

+ 2 (02) J 


Assembling Eqs. (82) — (85) in matrix form, we have 


' ' 

(fin) 

if (fin) 

1 

3 hp 
2 


T (io) 

fin) 

(?2 


1 

3 hp 
2 


. T (20) . 


' ' 
Q 3 

fin) 


1 

1 

1 

1 1 


1 1 

^ oT 
0 0 

fin) 


(82) 

(83) 

(84) 

(85) 

( 86 ) 

(87) 


Equations (86) and (87) relate the surface-averaged heat fluxes to the microvariables j. 

Next, we express the microvariables j* in terms of the surface-averaged temperatures in order 
to form the local conductivity matrix that relates the surface-averaged fluxes to the surface-averaged 
temperatures for a particular subcell (/3, 7). Hence, we proceed to define the surface-averaged tem- 
peratures as we have defined the surface-averaged heat fluxes above. At the left face of the subcell 
(/3, 7) i.e. at = —hp/2, the surface- averaged temperature T 2 is defined as 

T ~{P,1) = 1 f ll/2 T (A7) ( _^^(7) } ^(7) (88 ) 

Similarly, at the right face of the subcell (/?, 7), the surface- averaged temperature T 2 + ^’ 7 ^ is defined 
as 

t M0n) = j. /'’ /2 (89) 

h J-l-y / 2 Z 

At the bottom face of the subcell (/?, 7), the surface-averaged temperature T 3 ^ ,7 ^ is defined as 

tT (/?,7) = T~ r /2 T^\xf, J/)dx^ (90) 

h fi J-hf,/ 2 2 

At the top face of the subcell (/?, 7), the surface-averaged temperature is defined as 

A (/3,7) = A T^\xf\ l -^)dxf ] (91) 

h P J —hp/2 * 
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Substituting the expression for T given by Eq. (74) into the above definitions and performing 
the averaging procedure, we obtain 


rp (fill) 

1 2 

T".( 0 . 7 ) ^0 ^T.( 0 . 7 ) 1 ^0 rp(/3, 7) 

- i (00) 2 ( 10 ) 4 ( 2 °) 

( 92 ) 

rp+(fi->i) 

1 2 

T".( 0 . 7 ) 1 ^0 7) _J_ ^0 rp(/3, 7) 

- >o) + 2 J (10) + 4 J (20) 

( 93 ) 

rh~(fii 7) 

1 3 

71(0,7) 77(0,7) , 77(0,7) 

- >0) 2 V) + 4 (°2) 

( 94 ) 


7(0,7) 1 7 7(0,7) 1 7 7(0,7) 

- >0) + 2 V) + 4 (°2) 

( 95 ) 


Assembling Eqs. (92) — (95) in matrix form, we have 


r t + 1 

1 2 

(0,7) 

L A J 


r t + 1 

(0,7) 

LA \ 





T, 


( 10 ) 


T, 


( 20 ) 


(0,7) 


+ 


-too) 


T, 


Ti 


( 01 ) 
r ( 02 ) J 


1 (0,7) 


+ 


( 00 ) 


-( 00 ) 


T, 


(00) J 


(0,7) 


(0,7) 


Adding and subtracting Eq. (92) and (93), we obtain 


T (io) 

(0,7) 

" X 

7 

T 3 


' A ' 

(0,7) 4 

0 

. T (20) . 


1 

1 


La J 


. T (oo) . 


(96) 

(97) 


(98) 


Similarly, adding and subtracting Eq. (94) and (95), we obtain 


( 01 ) 

P (02) 


(0.7) 


J_ _ J_ 

\ \ 

l ^ 1 2 


A 

TT 


1 (0.7) 


12 

s 


0 

T (oo) 


(0.7) 


(99) 


Equations (98) and (99) relate the first and second order microvariables to the surface- aver aged 
temperatures and the zeroth order microvariables. In order to express the first and second order 
microvariables solely in terms of the surface- aver aged temperatures, we make use of the volume- 
averaged steady-state heat conduction equation. 


3.2.1 Heat conduction equation 


The heat flux in the material occupying the subcell (/?, 7 ) in the region 


L(/k7)| 


must satisfy the steady-state heat conduction equation given by Eq. (3) 
convenience, in a volumetric sense 


x 3 0,l) ^ 3*7> 
reproduced here for 


— 2^7h 


dq < 


(&7) 


dx . 




dq^P^ 

+ ^7-W =0 


dx 


*(7) 


(100) 


Substituting Eqs. (76) and (77) for the heat flux in the X 2 and x$ directions, respectively in the 
above equation, performing the required volumetric averaging and simplifying, we get 


[k 2 T m + k 3 T ( o2)] (/3 ’ 7) = 0 


( 101 ) 
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Substituting the second-order microvariables 
fying, we obtain 


and T, 


"( 02 )^ f rom Eqs. (98) and (99) and simpli- 


fy, 7 ) 


r (y,7) 

(00) 2 fc03.7) 


k^’^h 2 

(VO I i fD _ 


(A + T 2 -)(^) + A-^(T 3 + + T 3 -)(^> 


where 


2U^)l 2 


hi 


m 7 ) = fcf - 7) + 7) 

7 


Substituting Eq. (102) into Eqs. (98) and (99) and simplifying, we obtain 

~t~ i "1 09,7) 
r (io) 


- ( 20 ) 


[ (01) 
r (02) 


I 03,7) 



r t + i 
1 2 

03 , 7 ) 2k^'^ 

O 

O 

i 

r a ‘ 

L A J 

2 

[i iJ 

[a J 

A ‘ 

09 , 7 ) 2 A;^’ 7 ^ 

o 

o 

~ a ‘ 

L A \ 

k(^l 2 

1 1 

La J 


03,7) 


09,7) 


(102) 

(103) 

(104) 

(105) 


Equations (102), (104) and (105) express the microvariables explicitly in terms of the surface- 

averaged temperatures. Substituting Eqs. (104) and (105) into Eqs. (86) and (87) and simplifying, 
we obtain the local conductivity matrix as shown below 


where 


and 


Qt ' 

08,7) 

" Ml 

M2 

M3 

M4 

09,7) 

r T"+ i 
1 2 

(8,7) 

—Q 2 


Ml 

M2 

M3 

M4 


A 


Qt 


Ml 

M2 

M3 

M4 


T + 


. — $3 . 


_ Ml 

M2 

M3 

ft44 


La J 



(P> 7) _ *.(&7) _j _£(&7) 


m 


= ft, 


22 


„(0>7) _ 

A2 


,(Av) _ 

v 13 


_ J&7) 
Ml — M 


1 3 hpk^ r 

hp + 

3hpk:f r< ‘ > \ 
|2fc(y,7) ) 


.09,7) _ J/3,7) _ ,.03,7) _ ^hpk^^k^_ 


A4 


^23 


ft, 


24 


/2fc(A 7) 


ft 


(Av) 

33 


ft 


ft 


31 


= ft 


(&7) 

'44 


= -fe 


(Av) 


— + 


3fc, 


(0>7) 


Z 7 fc(^) 


(&7) _ 

12 


Ml — M 


3fc, 


(Av) 






(/5,7) _ ^(^,7) _ JAv) _ JAv) _ ^2 ’ ^ 


= ft. 


32 


= ft 


'41 


= ft 


42 


l yk^) 


(106) 


Therefore, we have formulated the local conductivity matrix which relates the surface- averaged 
heat fluxes to surface- averaged temperatures. The next step is to assemble the global conductivity 
matrix using the interfacial continuity conditions and the boundary conditions. 


NASA/CR— 2002-21 1909 


21 



3.2.2 Temperature continuity conditions 


The temperature continuity at the interfaces between adjacent subcells is applied in an average 
sense. Considering the j3 th interface, which is the interface between the subcells (/?, 7 ) and (/3+1, 7 ), 
the surface-averaged temperatures in the x 2 direction, T 2 + ^’ 7 ^ and T 2 must be equal. Hence, 

we can represent them using just one variable, i.e. 

f 2 + (/? ’ 7) = f 2 _(/3+1 ’ 7) = T 2 (/3+1,7) (107) 

Similarly, considering the interface, which is the interface between the subcells (/?, 7 ) and 
(/?, 7 + 1 ), and applying the temperature continuity in the x% direction 

f 3 +(/3,7) = t“ (/3 ’ 7+1) = T 3 (/3,7+1) (108) 

Equations (107) and (108) are similar to the temperature continuity conditions (19) — ( 22 ) in Section 
2 . Equations (107) and (108) hold true for (3 = 1, ..., Np — 1 and 7=1, ..., AT 7 — 1 , respectively. 
This gives rise to (Np — 1)TV 7 + (TV 7 — 1)7V# unknown surface-averaged temperatures defined at the 
subcell interfaces (both in X 2 and xs direction). The quantities 

rnO-1 7 ) /T : i(^Y# _ tb7) 

1 2 % 1 2 » ^3 5 3 

define the surface-averaged temperatures at the external boundaries of the composite. These quan- 
tities are either known or unknown depending on whether the temperature or heat flux is defined 
at the external boundaries. 


3.2.3 Heat flux continuity conditions 

The heat flux continuity at the interfaces between adjacent subcells is applied in an average sense. 
Considering the j3 th interface and applying the heat flux continuity in the x 2 direction 

Q+ (A7) - 0 2 (/3+1 ’ 7) = 0 (109) 

where Q^^’ 7 ^ is the surface-averaged heat flux going out of the right face (W 2 = hp/ 2) of the subcell 
(/?, 7) and Q 2 (/3+1,7) is the surface-averaged heat flux entering into the left face (X 2 = —hp+i/2) 
of the subcell (/? + 1 , 7 ) defined by Eqs. (79) and (78) respectively. Similarly, considering the ^y th 
interface and applying the heat flux continuity in the x% direction 


Q+(W _ Q-O^+i) = o 


(110) 


Equations (109) and (110) are similar to the heat flux continuity conditions (11) — (14). Using the 
local conductivity matrix (106), Eqs. (109) and (110) can be expressed in terms of the surface- 
averaged temperatures 


=H-(/3»7) 


, K (Prt)rf,-(Prt) 
“I" ^12 1 2 


' ™13 1 3 


1 1 J/3+b7)^+(/ ? +b7) , 

“T ^14 ^ 3 “T A € 21 1 2 t- 


(/3+l,7)^ — (/^H-1,7) 1 (^+1>7)7-i+(^+1»7) i (/3+l>7) n 

rboo -L 9 “I - 1190 -L q — I - Ai/q 4 -L q — U 


)rr~W,y) 1 ^(/?, 7)^+0, 7 ) 


(/3.7)t=--(/3,7) 
^34 J 3 


.(/3.7+l)T=i+0.7+l) 


k .(/?)7+1)t-> — (/5>7+l) 1 k .(/?)7+1)t->+(/3,7+1) i ^(/3>7+1)t->— (/?> 7+l) n 

ft , 42 -t 2 ft , 43 -t 3 ~T ^44 -£3 — u 


( 112 ) 
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Using the temperature continuity conditions given by Eqs. (107) and (108), Eqs. (Ill) and (112) 
can be simplified and written in terms of the common interfacial surface-averaged temperatures 


K 


(0,7) J,(0,7) 


12 

.(0,7 )r 

ns j 


+ (re- 


(0,7) 

li 


, (0+ 1,7) 

7 ^ ^22 


)Tr 


(0+1.7) k (0+1,7)^(0+2,7) 


21 


To 


, k ,(0.7)t=.(0.7) 
7” '*'14 


+ 


(0,7+1) _|_ k (0,7)t’(0+1,7) _l ^.(0+1,7) t-(0+1, 7+1) 


24 


'T 


+ ft, 


23 


= 0 


(113) 


K 


(0, 7)^(0, 7) _|_ k (0,7)^(0+1,7) _l ^ (0,7+1) t = >( 0, 7+1) _j_ t .(0,7+l)^r(0+l,7+l) 


32 
(0,7) ^ 


V 31 


+ «, 


42 


To' 


+ «. 


41 


+ 


K 34 ' ^3 


(0,7) 


+ 


l ,.(0,7) I -.(0,7+1) \y(0,7+l) I k .(0,7+1)t-'( 0>7+ 2 ) _ n 

1 1^33 \ ^44 / ^ 3 “T 1 + 4 q -L q — o 


M4 


M3 


(114) 


Thus, Eqs. (113) and (114) provide us with a total of (Np — l)iV 7 + ( N . 7 — 1)7V^ equations in terms 
of the common interfacial surface- aver aged temperatures and the surface- averaged temperatures at 
the external boundaries. 


3.2.4 Boundary conditions 

At the external boundaries of the composite, we have 2 (Np + AT 7 ) faces of the subcells where either 
the heat flux or the temperatures are defined. This gives rise to additional 2 (Np + TV 7 ) unknown 
surface-averaged quantities. The additional 2 (Np + N 7 ) equations are obtained from the imposed 
boundary conditions given by 

If' 7 ’ = T$,(x z ) (115) 

rf”' 71 = T^ ght (x 3) (116) 

where T^j t (x 3), and T^ ht (x 3) are the piece- wise uniform applied temperatures on the external 
boundaries in the X 2 — £3 plane. Similar reasoning holds for subcells (/?, 1), and (/?, iV 7 ) and the 
applied boundary conditions are given by 

tf 1} = 7ft L>2) (117) 

=T^{x 2 ) (118) 

If the heat flux at any of the boundaries is defined instead of temperatures, Eqs. (115) — (118) are 
modified accordingly. 

Note that, at least one and at least one at the external boundary need to be defined 
in order to prevent matrix singularity. If only heat flux is defined at the external boundaries, 
theoretically it gives rise to the possibility of having more than one temperature distribution. This 
is analogous to the rigid body motion if only tractions are defined at the external boundaries and 
not a single point is fixed in space. 

3.2.5 Solution for the surface-averaged temperatures 

Equations (113) and (114) together with the imposed boundary conditions (Eqs. (115) — (118)) 
provide us with the necessary (Np + 1)TV 7 + (TV 7 + 1 )Np relations for (Np + l)iV 7 + (N 1 + 1 )Np 
unknown surface- averaged variables (i.e., (Np — 1 )iV 7 + (N 1 — 1 )Np unknown common interfacial 
surface- averaged temperatures along with 2 (Np + N 7 ) unknown surface-averaged temperatures 
and/or heat fluxes at the external boundaries). 

The final system of equations is symbolically written as 

^ T = Q (119) 
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In the above equation, k is the global thermal conductivity matrix obtained after assembling 
the local conductivity matrices given by Eq. (106) using the local-global conductivity matrix 
approach as explained above. The matrix k essentially contains information on the geometry and 
thermal conductivities of the individual NpN^ subcells. The general format and assembly of the 
global thermal conductivity matrix k has been summarized in the next subsection. The vector 
T contains the unknown surface-averaged temperatures at the subcell interfaces and the surface- 
averaged temperatures at the outer edges of the composite (some of which are known), and is given 
by 

rp rrfi(l) rp(Nj) rp(l) rfiC^/s)] 

± ~ [■ l 2 > **•> ± 2 > x 3 > x 3 J 

where 

rp(v) _ ^(AT /3 + 1 , 7 ) j rp(/3) _ r fJiO-iP) + 1 ,/3) n 

2 L 2 ^ ^ 2 -I 3 L 3 ’ > 3 J 

The surface-averaged heat flux vector Q contains information about the piece-wise uniform heat 
fluxes defined at the external boundaries of the composite and consists mainly of zeros which are 
obtained after applying the interfacial heat flux continuity. It is given by 

Q = [Qy,..,Q^ ) ,Qy,...,Qf' 3) ] 


where 


= [Q^>,o,...,o,<3f» +1 ’ 7) ] 






Once Eq. (119) is solved for the surface- averaged temperatures at all the subcell interfaces and 
external boundaries, we substitute the surface- averaged temperatures back into Eqs. (102), (98) 
and (99), and obtain the microvariables which define the temperature field in each subcell. 


3.2.6 Assembly of the global thermal conductivity matrix k 

The general format and assembly of the global conductivity matrix k is summarized in this sub- 
section. k consists of four submatrices 

KXl K12 
^21 « 22 

where Kn and K 22 relate the quantities in their respective directions and have entries concentrated 
along the diagonal. The submatrices k\2 and K 21 represent the coupling of field variables in the X 2 
and x% direction and have entries scattered throughout. The global thermal conductivity matrix k 
is a square matrix of size [Np(Ny + 1) + Ny(Np + 1) x Np(N^ + 1) + + 1)]. The structure 

of the submatrix k n is shown below 



■ a^ i} 000000 0 

0 A y 00000 0 

0 0 .0000 0 

0 00.000 0 

K 11 = 

0 0 0 0.0 0 0 

0 0 0 0 0.0 0 

0 0 0 0 0 0. 0 

0 0 0 0 0 0 0 A^ 7) _ 
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where the size of k n is [N^(Np + 1) x N^(Np + 1)] . The structure of the submatrices is 

shown in the appendix. The structure of K 22 is similar to Kn. The structure of the coupling matrix 
K12 is shown below 


" b7 

b7 } 

°2 

Bf } 

g( 22 ) 

p> (2JV/3 ) 

°2 

-|-v (Ayi) 
^2 

pv (W 7 2) 

^2 

B (W-rJV/9) 


where the size of K12 is [N^(Np + 1) x Np(N^ + 1)] . The structure of the submatrices B^ 7 ^ is 
shown in the appendix. The structure of K 21 is similar to K12. 


3.3 Mechanical Reformulation 

The displacement field in the subcell (/?, 7) of the composite functionally graded in X 2 — £3 plane 
is approximated by the same second-order polynomial expansion in the local coordinates , and 
x^ as that given by Eqs. (32) — (34) in Section 2, reproduced below for convenience 


u 


(P,l) _ 


Xi£n 




_ w(^’7) 1 ^)i vd/^7) \ x^W^^ + -f3T^ 2 — jl _ L Jl\uAP^ 

U 2 — VV 2(00) ^ Xe 2 ^2(10) + ^2(01) + 2 ^ X 2 


-)^SS + 2 (3 "3 ^ 


iPa) 


(^’7) _L 7 ^( 7 ) (Pa) 


= W W)+* 


^3(10)+ *3 


3(01) 


+ §(3*f * - f )»3(2 0 ) 


hi 




+ - f )»' 3 '(02) 


Aha) 


For plane deformation (plane strain) case 


(120) 

( 121 ) 

( 122 ) 


£11 = 0 

and for generalized plane strain, e\\ is determined from the condition 


(123) 


N n Nt 


rhp /2 rl'y /2 
7 = 1 (3=1 J -hp/Z J—l'y / 2 


dx^ dx^ — 0 


(124) 


For an orthotropic elastic material occupying the subcell (/?, 7), the stress components are 
related to the strain components through the familiar Hooke’s law 


JAv) _ sy(P, i) (Prt) _ T(P,i) 

a ij ~ ° ijkl £ kl a ij 


(125) 


are the stiffness tensor elements for the material occupying the subcell (/?, 7), £\j’ v are the 
elastic strain components, and are the components of thermal stress vector given by 


Jfiri) 


7) _ -p(/^?7)^(/3,7) 

ij ij 


(126) 
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ut^) 

Fig. 5. Schematic of a subcell (/?, 7) showing the surface averaged tractions and displacements 

defined in the reformulated HOTFGM-2D. 


The local stiffness matrix relates the surface- averaged tractions to the surface- averaged displace- 
ments on the subcelhs faces as shown in Fig. 5. The surface- aver aged tractions are determined 
using the familiar expressions for tractions given in terms of the stress components 

M/3, 7 ) _ -(0,7) (J3,y) 

L i ~ a ji n 3 



where n^’ 7 ^ are the components of the unit vector normal to the face of the (/?, 7) subcell. At the 
left face of the subcell (/?, 7), the surface-averaged tractions t\ ^ ,7 ^ and 1 3 ^’ 7 ^ are defined as 


r2-C9,7) 

r 2 

r 3 


1 

hy 

1 



MPa) 

l 2 





blL jWwjW 

2 ) ax 3 

2 ’ x 3 ) ax 3 


(132) 

(133) 
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Similarly, surface-averaged tractions at the right face of the subcell (/?, 7 ) are defined as 


£ 2 t“ (ft 7 ) (^4 7 l v4 7) 

(134) 

111 2 “ 


f l "' 2 * 3 n(A 7 ) (^.4 7 l )<i4 7) 

(135) 

-£/2 2 



Surface averaged tractions at the bottom face of the subcell (/?, 7) are defined as 


x3— 0 , 7 ) 
l 2 

1 

hfiJ 

f hl 37/2 MPa)tAP) 
l 2 \ x 2 > 

-hp/2 


(136) 

r3-(/3, 7 ) 

r 3 

1 

hpJ 

f kl3 ^ 2 n (p) 

(3 V x 2 > 

-h.p/2 


(137) 


Surface averaged tractions at the top face of the subcell (/?, 7) are defined as 


f hel 2 


(138) 

-h p /2 



f’ 1 ' 2 


(139) 

~H! 2 

z 



Substituting Eqs. (131) , (125) , (126) , and (128) — (130) into Eqs. (132) — (139) , performing the 
required integration and assembling the resulting equations for the surface- averaged tractions in 
matrix form, we get 


j 2 + 

r 2 

p- 

ir, 


I (0.7) 


= c: 


(p, 7 ) 


22 


-j 3 / 1/3 

1 2 


W2(10) 

-j 3 /i/3 

2 


^2(20) 


i ( 0 . 7 ) 


+ eg- 7 ’ 


W-. 


3 ( 01 ) 


-w. 


(Ci20t\l + C , 22«22 + £ 23 ( 733 ) ^’ 7 ' ) 


n 

-To 


1 (/ 3 » 7 ) 


I (0.7) 


3 ( 01 ) 


J2+ 

r 3 

z 3 


(Pn) 


_ (o(/5,7) 


44 


1 

1 2 

_1 

1 2 




3 ( 10 ) 


W, 


3 ( 20 ) 


( 0 , 7 ) 


1 p(/3.7) 
-I- L-44 




2 ( 01 ) 


-w 


(0,7) 


2 ( 01 ) 


i3+ 

r 3 

f3- 

r 3 


(0,7) 




33 


1 3^, 

1 2 

3 Z-y 

1 “ 2 “ 


^3(01) 

^3(02) 


(0,7) 


+ 


^ 2 ( 10 ) 
— ^2(10) 


(£* 13(711 + £ 23<722 + £ 33 < 733 ) ^’ 7 ' ) 


-TT 


(0,7) 


£, 7 ) 


(140) 

(141) 


(142) 


f3+ 

l 2 

p- 

l 2 


I (Av) 


= a 


(Pn) 


44 


1 SLy 

, J 

1 2 




w, 


2 ( 01 ) 


n £. 7 ) 


2 ( 02 ) J 


1 ri(P,l) 
-T (- y 44 


Wi 


4£ 


3 ( 10 ) 


n £, 7 ) 


3 ( 10 ) J 


(143) 


Equations (140) — (143) relate the surface-averaged tractions to the first and second order mi- 
crovariables that determine the displacement field within the subcell (/?, 7). We need to relate these 
microvariables to the surface-averaged displacements, and surface- averaged temperatures obtained 
after solving the thermal problem, in order to formulate the local stiffness matrix. Hence, we 
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proceed to define the surface- averaged displacements. At the left face of the subcell (/?, 7), the 
surface- averaged displacements an d % are defined as 


- 2 - 03 , 7 ) _ 

1 

^ 7 /2 

.. 03 , 7 ) /_ 

h /3 

u 2 


— Z 7 /2 


2 

- _ 

i 

,/ 7 /2 

.. 03 , 7 ) 

h /3 

11 3 — 

^7 ^ 

— Z 7 /2 

u 3 l 

2 


Similarly, surface-averaged displacements at the right face of the subcell (/?, 7) are defined as 


: 2 + 03 , 7 ) _ 

1 

rl-y/2 

i • 

.. 03 , 7 )/ 

,hp 

~( 7 ) 

'2 


— / 7 /2 

ti 2 ( 

' 2 ’ 

^3 

: 2 + 03 , 7 ) _ 

1 

rl'y/ 2 

f • 

.. 03 , 7 )/ 

. h/3 


'3 


— / 7 /2 

tt 3 ( 

' 2 ’ 

^3 


Surface averaged displacements at the bottom face of the subcell (/?, 7) are defined as 


d“ (/3 ’ 7) = f / 

h 0 J-hp/2 2 

Surface averaged displacements at the top face of the subcell (/?, 7) are defined as 


j'hf 3/2 

4 ft,) ( 

> 0 ?) 
42 ’ 


(148) 

-^/2 





j-hp/2 

4 ft7l ( 

>09) 
, x 2 > 


(149) 

-ft/ 3/2 



z 



3+(/3,7) 

1 

r h p/ 2 

(43,7) / 

>09) 

7 

'2 

hp j 

/ 1 

—h/3/2 

« 2 ( 


7 

;3+03,7) _ 

1 

r h p/2 

..03,7)/ 

>09) 

7 

•3 

hp j 

/ 1 

f -^/ 2 

tt 3 ( 


7 


Substituting the expressions for given by Eqs. (121) and (122), and performing the above 

averaging procedure, we obtain 


,2-00,7) 


— — ^yp^’ 7 ) -u ^yp(^’ 7 ) 

^2(00) 2 ^ 2 ( 10 ) ^ 4 ^2(20) 


^2— (/3,7) _ ^T/T/0^’7) j 


3(00) 2 3 ( 10 ) 4 3(20) 


_2+(/3,7) _ j ^w(^»7) j ^TI/(^’7) 


2(00) ^ 2 2(10) ^ 4 fk 2(20) 


_2+(/3,7) jy(/?,7) _| — yE^’7) J ^.TI/(^’ 7 ) 


+ — w/ 


-3-(/3,7) 


3(00) ^ 2 3(10) ^ 4 '3(20) 

4 / 3 , 7 ) _ Ijwlftl ) 1 

2(00) 2 2 (01) 4 2(02) 


3-(/3, 7 ) _ w(/?,7) _ W(Al) , Sw(ft7) 

% - ' ' 3(00) 2 ”3(01) + 4 ''3(02) 

3+(/3, 7 ) _ ti/(/?,7) , 7 w(/5.7) , l J_ W ^) 

U 2 ~ ^2(00) + 2 ” 2 ( 01 ) + 4 ^2(02) 

3+(/3, 7 ) _ tiH0,7) 1 7) , 

U 3 ~ ^3(00) + 2 ”3(01) + 4 ^3(02) 
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Assembling Eqs. (152) — (159) for the surface- averaged displacements in matrix form 


u\ + 

,-.2- 


ul + 

- 2 - 


■Uo + 

d- 


(&7) 


( 0 , 7 ) 


( 0 , 7 ) 


2 

hp 
' 2 

h/3_ 

2 

hp_ 
' 2 


4 

4 J 


4 

4 


7 f 2 (io) 

( 0 . 7 ) 

4 - 

+ 2 ( 00 ) 

b / 2 ( 20 ) 

1 

Tl+oo) 


/J 


_k S 


4 J 


r .-73+ 1 (0,7) 




-3— 


L ^2 


^y ^ 

2 

"2 


^y /J 



W+io) 

( 0 , 7 ) 

+ 

W+oo) 


7E 3 (2o) 

7E 3 ( 00 ) 

" 

+1(01) 

( 0 . 7 ) 

+ 

7E 3 (oo) 1 


7E 3 ( 02 ) 

O 

O 


7E 2 (oi) 

( 0 . 7 ) 

+ 

5 

O 

O 


7E 2 (o2) 

O 

0 

1 


4 J L 


( 0 , 7 ) 


( 0 , 7 ) 


( 0 , 7 ) 


( 0 , 7 ) 


Adding and subtracting Eqs. (152) and (154), we get 


7E 2 (io) 

( 0 . 7 ) 

IE 2 (20) 



r J_ 

h i 

^/3 


J_ 1 

h l 

h< /3 



‘ t£+ ' 

( 0 , 7 ) 4 

0 1 


u\- J 


7 

0 

0 

1 


( 0 , 7 ) 


Adding and subtracting Eqs. (153) and (155), we get 


7E 3 (io) 

( 0 , 7 ) 

7E 3 (2o) 



r J_ 

jf 


J_ i 

h i 



®l + 1 

( 0 , 7 ) 4 

0 1 


[4- J 

h l 

7 

0 

0 

1 


( 0 , 7 ) 


Adding and subtracting Eqs. (157) and (159), we get 

1 ( 0 , 7 ) 


w. 


3(01) 


/ 2 

L i 


/ 2 


ul + 1 

4- 


W / 3(02) J 

Adding and subtracting Eqs. (156) and (158), we get 


( 0 , 7 ) 


72 

h 


0 

W+OO) J 


1 ( 0 , 7 ) 


IE 


2 ( 01 ) 

L ^2(02) J 


1 ( 0 , 7 ) 


1 1 


¥ 

S 


I 2 " 

r 


r ,.3+ i ( 0 , 7 ) 


L «2 


3- 


/2 

''7 L 


0 

^2(00) J 


4 ( 0 , 7 ) 


( 160 ) 

(161) 

(162) 

(163) 

(164) 

(165) 

(166) 
(167) 


Equations (164) — (167) relate the first and second order microvariables to the surface-averaged 
displacements and the zeroth order microvariables. In order to express the microvariables explicitly 
in terms of the surface-averaged quantities, we make use of the equilibrium equations. 


3.3.1 Equilibrium equations 


For the structure to remain in equilibrium, the stress field in the subcell (/3, 7 ) must satisfy the 
equilibrium equations. As in the original formulation of the higher-order theory, we satisfy the 
equilibrium equations, Eq. (36) in Section 2 , in an average sense, which is demonstrated by the 
following equations 


1 

+07) 




a^dx^dx^ 


= 0 , 


i = 2,3 


(168) 
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Substituting Eqs. (125) , (126) and (128) — (130) into Eq. (168) and performing the required 
integration yields, after simplification, the following equations 

3 (C 22 W 2m + C 44 W 2( oa))^ = (C 12 a n + C 22 a 22 + (169) 

3(^33^3(02) + C 44 W m 0 ) ) 6 o) = (Ci 3 an + 623022 + (170) 

Substituting the second order displacement microvariables from Eqs. (164) — (167) , we obtain 
the zeroth order microvariables in terms of the surface- averaged displacements and the first order 
temperature microvariables as shown below 


W 2?o’S = (^2 + + «2 ) (/3 ’ 7) + + U l ) (/3 ’ 7) - 7)T /fo7 


6,7) 


uirilfi-n) 


2 ( 00 ) 


2 a 


(fin) 


22 


012 n fin) 
7°22 


( 10 ) 


where 


and 


n fin) 12^ fin) 

<3 = +7< a i + + + +fey < a I + + a r) <A7) - a 33 ,) 41’ 7) 

ZUqq Z/T/^Oqq 


y 33 


')3'~ / 33 


a 


6 , 7 ) 


22 


= e 


6,7) , ^0 n ( 0 ’ 3 t ) 


22 


-h , 2 '-'44 

7 


/=((/?, 7 ) _ +/3,7) , ^7 +/3,7) 

°33 _ °33 "r ^2 0 44 


«22 


«33 = 


(( 7 x 20:11 + 6*22022 + 623033) ( ' /3,7 ' ) 


36* 


(0tl) 


22 


(613O11 + 623O22 + 633033)60) 


36 


6,7) 


33 


Substituting Eqs. (171) — (172) into Eqs. (164) — (167) and simplifying yields 


(171) 

(172) 

(173) 

(174) 

(175) 

(176) 


^2(10) 

6,7) 

1E 2 (20) 



2C, 


hfi 

{Pil) 


2 C\ 


h, 3 
(Pn) 


72 MPn) 
fc 7°22 

2Cfn ] \ 0 

12 r fin) 1 
S°22 l 


72 MPn) 
fc 7°22 


0 

1 



7:2 + 


u 2 


- 2 - 


L ^2 J 


O. 7 ) 


1 r ^3+ -1 6.7) 


u. 


-3- 


U. 


1 (&7) 


+ O 22 


-(10) J 


(177) 


^3(10) 

6,7) 

^3(20) 



1 


h(3 

9 

Z °33 

u2MP’l) 

+3° 33 


1 


hp 

9 /^<(Pi7) 
ZLy 33 

u2MPn) 

+3°33 



«3 + 1 


-2- 


% J 


6,7) 


o+/3,7) 
zu 33 
j. 2/6/3, 7 ) 

n 0^33 


O 

0 

1 


®r 1 

I' ij 


1 

1 

ICO CO 

153 


6,7) 


+ 


«22^ 


0 


1 6,7) 


T (oi) . 


(178) 
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^3(02) 


^°44 

h 2 MPrt) 

n (3^ 33 


u 2ri(Prt) 
n '/3 L '33 


^°44 

u2MPn) 

ri (3^ 33 

0 0 1 
1 1 


2+ -| (/3»7) 


112 ( 01 ) 


ZLy 22 


^22 

/2/^(^’3 / ) 

6 7°22 


o^C/^Tt) 

Z( ^22 


3+ 1 (At) 


-2+ -|(/3,7) 


Equations (177) — (180) relate the microvariables to the surface-averaged displacements and the 
first order temperature microvariables. Now we can obtain the local stiffness matrix for the subcell 
(/?, 7 ) by substituting the above equations into Eqs. (140) — (143) . Note that the terms involving 
the first order temperature microvariables can also be expressed in terms of the surface- averaged 
temperatures by making use of Eqs. (98) and (99). After simplification, we obtain 


t2+ 1 (0>7) 


ATn 

7712 

0 

0 

K is 

77i6 

K i 7 

77i 8 

7721 

7722 

0 

0 

7725 

7726 

K 27 

7728 

0 

0 

7733 

77 34 

7735 

77 36 

7737 

77 38 

0 

0 

7743 

77 44 

77 45 

77 46 

T7 47 

77 48 

i7 5 i 

7752 

7753 

77 54 

7755 

77 56 

0 

0 

77 61 

7762 

CO 

£ 

77 64 

77 65 

77 66 

0 

0 

T7 71 

T7 7 2 

T7 7 3 

T7 74 

0 

0 

T7 77 

T7 78 

77 8 i 

7782 

77 83 

77 84 

0 

0 

T7 87 

77 88 

r n 

r 12 

0 

0 1 

(At) 





r 21 r 22 0 

0 0 r 3 

0 0 r 4 

r 51 T52 0 

r6i r$2 0 

0 0 r 7 

0 0 r 8 


0 r 33 r 34 

0 r 43 r 44 
52 0 0 
62 0 0 

0 r 73 r 74 

0 r 83 r 84 


u 2 + 1 (/3 ’ 7) 


where 


tXP, 7 ) _ C 22 ’ 7 ' ) (A O ^ 22 ’ 

22 “^T ( d? 


-(2-3- 


iAAt) _ Ty-(P’l) _ 7AA7) 

21 16 — -“-25 — 21 26 


/7(/3,7) /2 

'~ / 22 S 


-<’ 7) = -4 ?’ 11 = 
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r4/5,7) _ 
-“■44 _ 


-( 4 - 3 - 


-( 2 - 3 ^ 


— _ — IsiPil) 

— ^45 — ^46 


jy'iPi'Y) _ 
^38 ~ 


isifli'y) _ 

iv 52 — 


js(Pn) 


AP, 7) _ &-(&7) _ 




w 23 "7 

3C 22 c^’ 7) ; 7 

/4/5,7) l2 
0 32 '7 


- _ &4& 7 ) — &4A7) 

— “63 — “64 


-( 4 - 3 - 


_ z>4/5,7) - 

— iV f^ — 


-(2-3^) 
v /40,7) ' 

Lyoo 


= -A 


7>4/5>7) _ 
“74 _ 


4/5,7) _ t^(/5,7) _ 


3C 33 c^’ 7) ; 7 

MPsi)u2 

°33 n (3 


ts(PA) 
^88 _ 


■(4-3^-r) 


_ r -(/ 5 , 7 ) _ 


^4/5,7) ^(4,7) 

733_/o _ o °33 

L [ 6 r^y 


4/5,7) _ o -4/3,7) / ^22 ’ 7 ' 


_ i-( 5.' ) _ 'J Q; 22 

1 1 o — 


3alo’ 7) C'oo’ 7) 


_-p(^>7) -p(^>7) _p(/V)0 

A 0/1 — A /IQ -L /I/I 


oAPAn^PA) 


_APA p(^>7) -n(A7) 

1 52 — 1 61 — 1 62 


-rg’ 7) = sag^(S|4 _ c<f>) 

r (A7) sag "'eg*' 

-L T'/i — 
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3.3.2 Displacement continuity conditions 


The displacement continuity at the interfaces between adjacent subcells is applied in an average 
sense. Considering the / 3 th interface, Fig. 6, the surface-averaged displacements in the X 2 direction, 
^2+(/3,7) an( j -2 (/3+ 1 , 7 )^ mus ^ p e equal. Hence, we represent them using just one variable, i.e., 

= u\~ {p+la) = -S2 (/3+1,7) (182) 

Similarly, applying the continuity of surface-averaged displacements in the x$ direction at the /3 th 
interface 

« 3 +(/3 ’ 7) = u ~ : 2(/3+1>7) = U 3 (/3+1 ’ 7) (183) 

Considering the 7 ^' interface, the surface-averaged displacements in the X 3 direction, and 

-3 (j)+ 1 . 7 ) ^ mus ^ equal and therefore can be represented using just one variable, i.e., 

_|+(/3, 7 ) = _3-(/3,7+l) = ^3(/3,7+l) (184) 

Similarly, applying the continuity of surface-averaged displacements in the X 2 direction at the ^y th 

_3+ ( /3, 7) = _3-(/3,7+l) = ^3(/3 >7 +l) (185) 

Equations (182) — (185) are similar to the displacement continuity conditions, Eqs. (59) — (62), 

in Section 2 . Equations (182) — (183) and (184) — (185) hold true for /3 = 1, ...,Np — 1 and 7 = 

1 , ..., N 1 — 1 , respectively. This gives rise to 2(Np — 1)AT 7 + 2 (TV 7 — 1 )Np unknown surface- averaged 
displacements defined at the subcell interfaces (both in X 2 and X 3 directions). The quantities 


2(1,7) -2(1,7) 2(^+1, 7 ) 

a 2 5 a 3 > a 2 


2(^+1, 7 ) 

u 3 


, U. 


3(/3,l) -3(/3,l) 


3(AW 7 +1) 3(AW 7 +1) 

a 2 ? a 3 


define the surface-averaged displacements at the external boundaries of the composite. These 
quantities are either known or unknown depending on whether tractions or displacements are 
defined at the external boundaries. 


ji+ifrr) + j^2-(j0+i ,r) _ Q = ^-(je+Lr) _ 



Fig. 6. Traction and displacement continuity applied at the interface of the subcells in an average 

sense. 
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3.3.3 Traction continuity conditions 

The traction continuity at the interfaces between adjacent subcells is applied in an average sense. 
Considering the j3 th interface, Fig. 6, the continuity of tractions in the X 2 direction is ensured by 

? 2+(/3, 7 ) + f 2-(/?+l >7 ) = q ( 186 ) 

where t is the surface-averaged traction at the right face ( 0 C 2 = hpj 2) of the subcell (/?, 7) and 
t 2 ^ +1,7 ^ is the surface-averaged traction at the left face (#2 = ~hp+ 1/2) of the subcell (/? + 1, 7) in 
the X 2 direction, defined by Eqs. (89) and (88), respectively. Similarly, the continuity of surface- 
averaged tractions at the f3 th interface in the X 3 direction is ensured by 

t -2+(/3 >7 ) + fM/3+1,7) = q (187) 

Considering the j th interface, and applying the continuity of tractions in the X 2 direction in an 
average sense we obtain 

t -3+(/3, 7 ) + 73-(/3,7+l) = q ( 188 ) 

The continuity of surface-averaged tractions in the x 3 direction at the r y th interface is ensured by 

^-3+0,7) + 73-0.7+1) = o ( 189 ) 

Equations (186) — (189) are similar to the traction continuity conditions, Eqs. (51) — (54), in 
Section 2. Using the local stiffness matrix (181), Eqs. (186)-(189) can be expressed in terms of the 
surface- aver aged displacements 

(K n ul + + k 12u\- + K 15 u 3+ + K 16 u 3 ~ + K 17 u 3+ + K 18 u 3 ~)^+ 

{K 21 u 2+ + K 22 u 2 ~ + K 25 u 3+ + K 26 u 3 2 ~ + K 27 u 3+ + K 28 u 3 -)^ +1 ^ = 

(r n f 2 + + r 12 f 2 -)^) + (r 2 if + + t 22 t 2 )^ +1 ^ (190) 

(i^ss^l + + K 34 u 2 ~ + K 3 b u 3+ + K 36 u 3 ~ + K 37 u 3+ + K 38 u 3 ~)^+ 

(K 43 u 2+ + K 44 u 2 ~ + K 4 b u 3+ + K m u 3 ~ + K 47 u 3+ + K 48 u 3 ~)^ +1 ^ = 

(T 33 T+ + r 3 4T 3 _ ) (/3 ’ 7) + ( r 43T 3 + + r 4 4r 3 _ ) (/3+1,7) (i9i) 

{K bl u 2+ + K b 2 u\~ + K b 3 u 2+ + K b 4 u\~ + K bb u 3+ + K b 6 u 3 ~)^+ 

( K ei ul + + K 62 u 2 ~ + K 63 u 2+ + K m u\~ + K 65 u 3+ + K 6 6 ^-) (/3 ’ 7+1) = 

(r 5 iT 2 + + r 52 T 2 _ ) (/3 ’ 7) + (r 6 i t+ + r 62 f 2 )^+v (192) 

(K n u 2+ + K 72 u\- + K 73 u 2+ + K 74 u 2 ~ + K 77 u 3+ + K 78 u\-)^+ 

(K 8 iu 2+ + K 82 u\~ + K 83 u 2+ + K 84 u\~ + K 87 u 3+ + K 88 u 3 ~)^ +1) = 

(r 73 f 3 + + r 74 f 3 -)^ + (r 83 f 3 + + r 84 t 3 “) (/3 ’ 7+1) (193) 

Using the displacement continuity conditions given by Eqs. (182) — (185), Eqs. (190) — (193) can 
be written in terms of the common surface- averaged displacements 

^(J 7) _20,7) + (i f0,7) + x O+l,7)^2(/3+l,7) + ^0+1,0-20+2,7) + ^0,7) ^0,7) + 

1) i T?(P+ 1 >'y)j-MP+ 1 >'y) , t^(/ 3+1,7) 9 -3(/3+1,7+1) , 
iV 15 a 2 “r iV 26 a 2 “r iV 25 a 2 ' ^18 a 3 ' 

^ (A7) ^3 (/3, 7+1) ^^+-*-’T)^3(/3+1,7) ^-^+1,7) -3(/3+l,7+l) _ 

r 0.7)^0, 7) + (r 0,7) + r 0+l,7) ):f 0+l,7) + rg+ 1 , 7 ) jO+2,0 (194) 
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K. 


0 3 ’T')„- 7 2(/3,7) 


34 


+ (K. 


( 0 , 7) 


33 


+ AT 


03 + 1 , 7 ) x - 209 + 1 , 7 ) 


44 


+ A' 


03 + 1 , 7 ) - 209 + 2 , 7 ) 


43 


M. 


+ A. 


(0^) r M0n) 


36 


M. 


a:: 


(/3+l>7)„-;3(/3+l,7) 


46 


U, 


+ K. 


(A7)„- 7 3(/3,7+l) 


35 


U. 


'(/3+l,7)^3(/3+l,7+l) 


1^03+1.7) .-j309+1, 7 ) 

11 48 % 


r 


(/ 3 , 7 ) 

34 


7 “>(/ 3 > 7 ) 1 p(/ 3 > 7 ) 
J 3 + 1 33 


7 i(/ 3 , 7 +l) 1 ■p 
J 3 + x 44 


+ -^45 “2 

+ K^ ] ul^ +l) 

(/ 3 +l, 7 ) 


+ at: 




38 


U 


+ 

+ 


+ a; 


09+1, 7)-303+l, 7+1) _ 


47 


■u. 


7 >(/ 3 +l, 7 ) 1 -p(/ 3 +l, 7 )T->(/ 3 +l, 7 +l) 
1 Z +J -43 1 Z 


( 195 ) 


K, 




52 


A, 


(A 7 )„-; 2 (/ 3 , 7 ) 


54 


+ K, 

+ A, 


(^’7)»-7 2 (/3+!.7) 


51 


U. 


(/3>7)r7 2 (/3+l,7) 


53 


It. 


a:, 


(/ 3 >7)^3(/3,7) 


56 


U 


+ AT, 

+ A', 

03 > 7 ) 


(/ 3 > 7 +l)„- 7 2 (/ 3 , 7 +l) 
62 “2 

(/ 3 , 7 +l)- 2 (/ 3 , 7 +l) 
64 “3 


+ at ( 


( y 9 , 7 +l)_ 2 ( y 9 +l, 7 +l) 


61 


U. 


+ ^>+ Ki 


+ A, 

K/ 3 > 7 +l)\- 30 ?, 7 +l) 


M. 


66 




09 , 7 +l)_ 203 +l, 7 +l) 
63 

+ AT, 


+ 

+ 


(/ 3 , 7 +l).- 7 3 (/ 3 , 7 + 2 ) _ 


65 


U 


0.7) jiO.7) _j_ p0.7)^r0+1.7) _j_ p0.7+1)t>0.7+1) _l pO, 7+1)7^03+1, 7+1) 


52 


+ r 


51 


'T 


+ r, 


62 


To 


+ r, 


61 


'To 


(196) 


h TPn)-XPa) 

tv a % 


A$’ 7) u^’ 7) + A'g’ 7) u^ (/3+1 ’ 7) + A^’ 7+1) ^ (/3 ’ 7+1) 

3 , 7 ) _i_ 1 ^ 03 , 7 ) ,7203+1,7) , jpO 3 ’ 7 +l), 7 2 (/ 3 , 7 +l) 
+ -fv 73 U 3 + -t7 8 4 % 

K (P,l) j-MPrt) , (TsiP’T 4 . / p (/ 3 > 7 + 1 )' ) 1 3 (/ 3 / 


+ at^’ 7+1) ^ 




(197) 


Thus, Eqs. (194) — (197) provide us with a total of 2 (Np — 1)TV 7 + 2 (TV 7 — 1)7V# equations in terms 
of the common interfacial surface- averaged displacements and the surface- averaged displacements 
at the external boundaries. 


3.3.4 Boundary conditions 

At the external boundaries of the composite, we have 2 (Np + TV 7 ) faces of the subcells where either 
tractions or displacements are defined. This gives rise to additional 4 (Np + AT 7 ) unknown surface- 
averaged quantities. The additional 4 (Np + N^) equations are obtained from the imposed boundary 
conditions given by 


(198) 

‘P , ’ ) = ‘Su(*3l 1 = 2.3 (199) 

where tiiftto ) and t right( x 3 ) are piece- wise uniform surface tractions applied on the vertical bound- 
aries in the X 2 — £3 plane. Similar reasoning holds for subcells (/?, 1 ), and (/?, TV 7 ). Alternatively, 
if the displacements are specified on the vertical boundaries, then the applied boundary conditions 
are given by 


u z(/3,i) = u (P) tom (x 2 ) (200) 

^ (W = 4S(® 2 ) * = 2 ,3 (201) 

For other type of boundary conditions, the Eqs. (198) — (201) are modified accordingly. 

Note that at least one and at least one should be defined at the external boundary in 
order to prevent rigid body motion. 
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3.3.5 Solution for the surface-averaged displacements 

Equations (194) — (197) together with the imposed boundary conditions (198) — (201) provide us 
with the necessary 2(Np + l)N^ + 2(Ny + l)Np relations for the 2(Np + l)Nn f + 2(N^ + l)Np unknown 
surface-averaged variables, i.e., 2(Np — l)iV 7 + 2(AT 7 — 1 )Np unknown common interfacial surface- 
averaged displacements along with 4 (Np + TV 7 ) unknown surface- averaged displacements and/or 
surface tractions at the external boundaries. Also, the unknown uniform strain En is determined 
from the generalized plane strain condition given by Eq. (124). For plane strain, En is zero. 

The final system of equations obtained is symbolically written as 

K U=t (202) 


In the above equation, K is the global stiffness matrix obtained after assembling the local stiffness 
matrices given by Eq. (181) using the local-global stiffness matrix approach as explained above. 
The matrix K essentially contains information on the geometry and thermomechanical properties 
of the individual NpN 1 subcells. The general format and assembly of the global stiffness matrix K 
has been summarized in the next subsection. The vector U contains the unknown surface- averaged 
displacements at the subcell interfaces and the outer edges of the composite and is given by 


where 


U= [uft... 

r 2(JV 7 ) 

2(1) 2 (N 7 ) 

u 3 ; -? u 3 5 

3(1) 3 (N P ) 3(1) 

u 2 > u 2 ) u 3 

M) = 1 

, ..., u 3 ,£llj 

«l w = 

[a 

2(^+1, 7 )i 

>•••>“ 2 J 


2(^+1, 7 ) ] 
u 3 \ 

^ = 


3(^+1, /3) 1 

? "2 J 

4 w = i4 {1 ’ 

-3(N/3+l,{3)-\ 
> u 3 \ 


The surface-averaged traction vector t contains information on the applied boundary conditions 
and the mechanical effects produced by thermal loading. 

Once Eq. (202) is solved for the surface-averaged displacements at all the subcell interfaces and 
external boundaries, we substitute the surface- aver aged displacements back into Eqs. (171) — (172) 
and Eqs. (164) — (167) and obtain the microvariables W^^n) which define the displacement field 
in each subcell. 


3.3.6 Assembly of the global stiffness matrix K 

The general format and assembly of the global stiffness matrix K is summarized in this subsection. 
K consists of eight submatrices 



Kn 

0 

K 13 

Ku 

0 

k 22 

k 23 

k 24 

K 31 

K 32 

K 33 

0 

K 41 

K 42 

0 

k 44 


where Kn,K 22 ,K 33 and K44 relate the quantities in their respective directions and have entries 
concentrated along the diagonal. The remaining submatrices represent coupling of the field variables 
in the X 2 and X 3 directions and have entries scattered throughout. The size of the global stiffness 
matrix K is [2 Np(N 7 + 1) + 2 N 7 (Np + 1)] x [2Np(N 7 + 1) + 2 N 7 (Np + 1)] for plane strain. An 
additional row and column is added for the generalized plane strain case. The structure of the 
submatrix Kn is shown below 
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0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

. 0 0 

0 . 0 

0 0 A^ 7) . 

where the size of Kn is [Ny(Np + 1) x Nj(Np + 1 )] . The structure of the submatrices is 

shown in the appendix. The structure of K22, K33 and K44 is similar to Kn. The structure of the 
coupling matrix K13 is shown below 


B( n ) 

°13 

b( 12 ) 

**13 

p>(lAW 1 

b 13 

b( 21 ) 

a 13 

b( 22 ) 

**13 

13 13 

■ 0 CW 7 I) 

_ **13 

-d(N 7 2 ) 

^13 

B 13 


Kn = 


A (1) 

A n 

0 

0 

0 

0 

0 

0 

0 


v ( 2 ) 

Mi 

0 

0 

0 

0 

0 

0 


0 0 0 

0 0 0 

. 0 0 

0 . 0 

0 0 . 

0 0 0 

0 0 0 

0 0 0 


where the size of K13 is [ A 7 ( Np + 1) x Np(N 1 + 1)] . The structure of the submatrices is 

shown in the appendix. The structure of the remaining coupling matrices is similar to K13. 

4 Mesh-Sensitivity and Validation Studies 

The two-dimensional formulation of the original higher-order theory was discussed briefly in Section 
2 . In Section 3 , an efficient reformulation of the higher-order theory was developed and discussed 
in detail. The procedure for determining the various field quantities using the reformulated higher- 
order theory was also outlined in Section 3 . The next step is to verify the efficiency and accuracy 
of this reformulated version for various thermal, mechanical, and thermomechanical problems. 
In this section, we investigate the convergence of thermal and mechanical field quantities with 
mesh refinement and also verify their accuracy upon comparison with analytical and finite-element 
solutions. 

4.1 Mesh Sensitivity: Thermal Problem 

In this subsection, the convergence of temperature field with mesh refinement is investigated and 
the results are compared with an analytical solution. The problem definition, investigated geometry 
and meshing are shown in Fig. 7 . As shown in Fig. 7 (a), the cross-section of the block has unit 
dimensions in the re 2 — a; 3 plane. The dimension of the block is considered infinite in the out-of-plane 
(aq) direction which, however, does not play a role in the temperature field analysis. The block is 
subjected to a temperature of 100°C', which is held constant, at the left, top and bottom faces and 
a temperature of 200 °C at the right face. As shown in Figs. 7 (b), (c) and (d), the cross-section of 
the block in the aq — x ?> plane is discretized into 4 x 4, 12 x 12, and 32 x 32 subcells, respectively, 
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Problem Definition 


•*3 Temperature B.C.'s at the edges 
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Homogeneous Material 


k= 25 W/m-°C 
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12x12 subcells 
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Fig. 7. Problem definition, investigated geometry and discretization used for mesh sensitivity 

studies: thermal problem. 
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of uniform size. The material within the block is homogeneous and the thermal conductivity is 
taken to be 25 W/m — °C. However, the analytical solution of the Laplace’s equation indicates that 
the temperature field is independent of the thermal conductivity for the homogeneous case. 

For the given boundary conditions, geometry and meshing, the temperature field was generated 
using the reformulated higher-order theory. Analytical solution for the temperature field was ob- 
tained by solving the Laplace equation using the standard Fourier series approach, Zhong (2002). 
The contour plots of the temperature fields for the 4 x 4 , 12 x 12, and 32 x 32 subcell meshes are 
shown in Figs. 8 (a), (b) and (c), respectively, and are compared with the analytical solution shown 
in Fig. 8 (d). As observed from the contour plots, the 4x4 subcell mesh approximates the actual 
temperature field only in a rough sense. As the mesh is refined to 12 x 12 subcells, the temperature 
field becomes almost identical with the actual temperature field except near the ( 1 , 0 ) and ( 1 , 1 ) 
coordinates in the x 2 — X 3 plane. These are the points of temperature discontinuity and hence 
greater mesh refinement is required in order to properly capture the actual temperature near these 
points. This is achieved with the 32 x 32 subcell mesh for which the temperature field, Fig. 8 (c), is 
visually identical to the actual temperature held, Fig. 8 (d), obtained from the analytical solution. 

Next, we consider the convergence behavior along several cross-sections. For the 4x4 subcell 
mesh shown in Fig. 7(b), the cross-section along the line X 3 = 0.25 happens to be the interface 
between the subcells ( 1 , 7 ) and ( 2 , 7 ), and the cross-section along the line £3 = 0.5 is the interface 
between the subcells ( 2 , 7 ) and ( 3 , 7 ). Similarly, for the 12 x 12 and 32 x 32 subcell cases, the 
cross-sections along the lines X 3 = 0.25 and £3 = 0.5 run along the corresponding subcell interfaces. 
The pointwise temperature distribution (also displacements and stresses) along each interface can 
be calculated using the microvariables associated with the subcell on either side of the interface. In 
general, the microvariables belonging to subcells on the opposite sides of the interface are different. 
Therefore, the pointwise interfacial temperatures calculated using the microvariables on either side 
will be different since the thermal/heat flux continuity conditions across the interface are applied 
in a surface-average sense. However, with the refinement in mesh, the interfacial temperatures 
calculated using the subcell microvariables on either side of the interface should converge. The 
exception occurs when the interfacial line happens to be the line of symmetry. In that case, the 
magnitude of the subcell microvariables on the opposite sides of the interface is identical and hence 
the interfacial temperatures. 

Figure 9 shows the temperature distributions along the lines £3 = 0.25 and X 3 = 0.5. The 
symbol X 2 denotes the common interfacial surface-averaged temperatures which are the basis of 
the reformulation. The product of the global thermal conductivity with these interfacial surface- 
averaged temperatures yields the applied thermal/heat flux boundary conditions. T + denotes the 
interfacial temperature obtained using the microvariables of the subcells lying below the interfacial 
lines. T~ denotes the interfacial temperature obtained using the microvariables of the subcells 
lying above the interfacial lines. As observed in Fig. 9 (a), the temperature distributions along 
the line X 3 = 0.25 calculated using the subcell microvariables on the opposite sides of the interface 
are quite different for the 4x4 subcell case. However, the interfacial temperature distributions 
tend to converge to the same values for the 12 x 12 subcell case and practically coincide for the 
32 x 32 subcell case, Figs. 9 (c) and (e), respectively. Also, the surface- averaged temperatures (T 2 ) 
tend to converge in a piece- wise uniform manner to the actual temperature distribution with mesh 
refinement. For the cross-section along the line X 3 = 0.5, the temperature distributions, T + and 
T _ , calculated using the subcell microvariables on the opposite sides of the interface are identical 
even for the 4x4 subcell case, Fig. 9 (b). This is because the line X 3 = 0.5 happens to be the 
line of symmetry for the given boundary conditions and hence the microvariables belonging to the 
subcells on the opposite side of the interface are identical. This provides an additional verification 
of the results obtained using the reformulated higher-order theory. 
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Fig. 8. Temperature field for mesh discretizations and boundary conditions given in Fig. 7 and 

comparison with analytical solution. 
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Fig. 9. Convergence studies for temperature distributions along the cross-sections X 3 = 0.25, 

shown in Fig. 7. 


NASA/CR— 2002-21 1909 


41 


Figure 10 shows the temperature distributions along the lines x 2 == 0.5 and X 2 = 0.75. These 
cross-sections do not lie along the line of symmetry and hence the corresponding interfacial distribu- 
tions exhibit trends similar to those observed in the cross-section along the line x% = 0.25 described 
above. The temperature distributions T + and T~ calculated using the subcell microvariables on 
the left and right sides of the interface, respectively, are quite different for the 4x4 subcell case, 
Figs. 10 (a) and (b). These interfacial temperature distributions tend to converge to the same 
values for the 12 x 12 subcell case, Figs. 10 (c) and (d), and coincide for the 32 x 32 subcell case, 
Figs. 10 (e) and (f). Also, the surface- aver aged temperatures (T 3 ) tend to converge to the actual 
temperature distribution with mesh refinement. 

Figure 11 shows the temperature distributions along the boundaries X 3 = 0 and X 2 = 1. As 
observed in the figure, these distributions calculated using the reformulated higher-order theory 
approach the applied boundary conditions with mesh refinement. However, they do not exactly 
match the applied boundary conditions due to the presence of the temperature discontinuity at 
the (1,0) and (1,1) coordinates. Virtually, an infinitely dense mesh near these points would be 
required to exactly match the applied boundary conditions. Even in the case of the analytical 
solution, a large number of terms in the Fourier series expansion is needed in order to obtain 
converged solution near the points of temperature discontinuity. At the point of discontinuity, 
however, only the average value is obtained according to the well-known theorem. 

The interfacial temperature distributions calculated using the subcell microvariables on the 
opposite sides of the interface and averaged at each point across the interface are plotted along the 
various cross-sections (along the lines X 3 — 0.25, 0.5 and X 2 — 0.5, 0.75) in Fig. 12. These plots were 
generated for the different meshes considered above and the results are compared with the analytical 
solution. As observed in these cross-sectional plots, averaging the temperatures calculated using 
the subcell microvariables on the opposite sides of the interface produces acceptable results even for 
the rough mesh as in the 4x4 subcell case considered here. The difference between the analytical 
solution and the 4x4 subcell mesh is greater in the regions of high temperature gradients. Since the 
considered cross-sections are removed from the points of temperature discontinuity, the temperature 
distributions generated using 12 x 12 and 32x32 subcell meshes coincide with the analytical solution. 

4.2 Validation: Thermal Problem 

In the case of functionally graded materials, the microstructural gradation is typically varied gradu- 
ally in a manner that depends on the boundary conditions in order to obtain the required optimized 
composition profile. This results in continuous or discrete gradation of microstructure as described 
in Section 1 , and therefore continuously or discretely varying properties. In this subsection, a 
heterogeneous composition with discrete variation of microstructure is considered. In order to 
generate discretely varying thermal conductivity fc, a continuous (exponential) function of spatial 
coordinates is employed as a basis. The thermal conductivity is then calculated at the center of 
each subcell according to the given function and is assumed to be constant within the subcell. 
The results are compared with the finite-element solution obtained using ANSYS, and the effect of 
thermal conductivity variation on temperature field and its gradients is discussed. 

Four different thermal conductivity variations were considered. The problem definition, inves- 
tigated geometry and meshing are given in Fig. 13. The investigated geometry and meshing are 
the same for all four cases. The cross-section of the block has unit dimensions in the X 2 — X 3 plane 
and is considered infinitely long in the out-of-plane (aq) direction. As observed in Subsection 4.1, 
the temperature field for the homogeneous case converged fairly well with the 12 x 12 subcell mesh 
except at the points of temperature discontinuity, and it was identical to the actual solution for 
the 32 x 32 subcell mesh. Therefore, in the cases considered here, the cross-section was discretized 
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Fig. 10. Convergence studies for temperature distributions along the cross-sections X 2 = 0.5, 0.75 

shown in Fig. 7. 
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Fig. 11 . Convergence studies for temperature distributions along the edges X 3 = 0 and X 2 = 1 

shown in Fig. 7. 
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Fig. 12. Convergence studies for averaged interfacial temperature distributions along the 
cross-sections shown in Fig. 7, and comparison with analytical solutions. 
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into 25 x 25 subcells of uniform size. In the first case, Fig. 13 (a), the block is subjected to a 
temperature of 100°C, which is held constant, at the left, top and bottom faces and a temperature 
of 200 °C at the right face. The thermal conductivity k ( W/m — °C ) is assumed to be exponentially 
increasing in the x 2 direction according to the function 

k(x 2 ) = ie 5 * 2 (203) 

In the second case, Fig. 13 (b), the applied temperature boundary conditions remain the same but 
the thermal conductivity k is assumed to be exponentially decreasing in the X 2 direction according 
to the function 

k(x 2 ) = 60e~ 5x2 (204) 

In the third and fourth cases, the block is subjected to a temperature of 100° ( 7 , which is held 
constant, at the left and top faces, and a temperature of 200°C at the right and bottom faces. 
The thermal conductivity k is assumed to be vary exponentially in both the X 2 and X 3 directions 
according to the functions 


k(x 2 ,* 3 ) = 10e 

(205) 

k(x 2 ,xs) = 10 e 2x2 ~ 2x3 

(206) 


4.2.1 Reformulated HOTFGM and finite-element comparison 

The temperature field calculated using the reformulated version of the higher-order theory is com- 
pared with the finite-element solution. The finite-element solution was obtained by simulating the 
above cases in ANSYS using 8-node thermal elements called Plane77. The mesh discretization used 
in the finite-element analysis is the same as that used in the reformulated higher-order theory. Also, 
the Plane77 element in ANSYS has the same order (quadratic) of temperature field approximation 
as in the reformulated higher-order theory. For the first case, Fig. 13 (a), the temperature field 
contour plots obtained using the reformulated higher-order theory and ANSYS are shown in Figs. 
14 (a) and (b), respectively. For the second case, Fig. 13 (b), the temperature field contour plots 
obtained using the reformulated higher-order theory and ANSYS are shown in Figs. 14 (c) and 
(d), respectively. For the third case, Fig. 13 (c), the temperature field contour plots obtained using 
the reformulated higher-order theory and ANSYS are shown in Figs. 15 (a) and (b), respectively. 
Finally, for the fourth case, Fig. 13 (d), the temperature field contour plots obtained using the 
reformulated higher-order theory and ANSYS are shown in Figs. 15 (c) and (d), respectively. 

As observed in the contour plots, the temperature fields obtained from the reformulated higher- 
order theory and the finite-element analysis match very closely. However, at the points with the 
coordinates (1,0) and (1,1), which are the points of temperature discontinuity, the temperature 
field obtained from reformulated higher-order theory shows better convergence than the finite- 
element solution. This is because the boundary conditions in the higher-order theory are applied in 
a surface-average sense, whereas in the finite-elements the boundary conditions are applied at the 
nodes. Therefore, higher mesh refinement is required at the points of temperature discontinuity in 
the finite- element case for this particular choice of element. 

4.2.2 Effect of thermal conductivity variation 

According to the Fourier’s law of heat conduction, the heat flux is in the direction of the nega- 
tive temperature gradient. Moreover, the temperature field is modulated by the variation in thermal 
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Fig. 13. Different thermal conductivity variations and mesh discretizations used for validation 

studies: thermal problem. 
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Fig. 14. Temperature field comparison obtained using HOTFGM and FEA for thermal 

conductivity variations in X 2 direction. 
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Fig. 15. Temperature field comparison obtained using HOTFGM and FEA for thermal 
conductivity variations in both X 2 and X 3 direction. 
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conductivity. Increasing the thermal conductivity in a particular direction shifts the temperature 
field towards the opposite direction, thereby decreasing the temperature gradient and vice-versa. 

The temperature field for the same geometry and boundary conditions as in the first and the 
second case discussed in this subsection, but for the homogeneous material, was presented in Fig. 8 
(d). In the first case, increasing the thermal conductivity in the X 2 direction shifts the temperature 
field in the negative X 2 direction, Fig. 14 (a), relative to the homogeneous case. In the second case, 
on the other hand, decreasing the thermal conductivity in the X 2 direction shifts the temperature 
field in the positive X 2 direction, thereby increasing the temperature gradient near the right face 
of the block, Fig. 14 (c). Similar behavior is observed for the two-dimensional variation of thermal 
conductivity shown in Figs. 15 (a) and (c). 

4.3 Validation: Mechanical Problem 

In Subsections 4.1 and 4.2, mesh sensitivity and validation studies were conducted for thermal 
cases with both homogeneous and heterogenous materials. In this subsection, the classical Eshelby 
problem for the pure mechanical loading is considered. This problem involves an elliptical fiber 
inclusion in an infinite matrix with uniform surface tractions applied over the boundaries at infinity. 

Here, for the purpose of comparison, a circular inclusion is considered instead of an elliptical 
one and plane strain analysis is carried out. The matrix is assumed to be made up of epoxy and 
its cross-section has unit dimensions in the X 2 — x% plane. A very small glass fiber is embedded in 
the matrix so that the matrix practically behaves as infinite and the effects of the fiber inclusion’s 
presence on stresses are not felt near the edges. The material properties of glass and epoxy used 
in the analysis are listed in Table 1. The matrix is subjected to normal surface tractions of unit 
magnitude in the X 2 direction at the outer faces as shown in Fig. 16 (a). 


Material 

E ( GPa ) 

V 

Glass fiber 

69.0 

0.2 

Epoxy matrix 

4.8 

0.34 


Table 1. Material properties of constituent fiber and matrix phases. 

Because of the type of boundary conditions and the investigated geometry, the problem is 
symmetric about the cross-sections along the lines X 2 = 0.5 and xs = 0.5. The mid-points along 
the lines X 2 = 0, 1 and xs = 0, 1 are not expected to move in the xs and X 2 directions, respectively. 
Therefore, in order to prevent rigid body motion, the middle two subcells along the lines X 2 = 0 
and xs = 0 were constrained from moving in the xs and X 2 direction, respectively. 

The full mesh used in the reformulated higher-order theory is shown in Fig. 16 (a). In order 
to capture the high stress gradients near the interface of the glass fiber and epoxy matrix, which 
occur due to the large material property mismatch between the two materials, very refined mesh 
was used in the interface’s vicinity. In order to clearly see the refined mesh near the inclusion, the 
magnified mesh discretization is shown in Fig. 16 (b). 

Stress contours obtained using the reformulated higher-order theory are shown in Figs. 17 
(a), (c) and (e) and compared with the analytical solution (cf., Dugdale and Ruiz (1971)) in 
Figs. 17 (b), (d) and (f), respectively. As observed from the contour plots, the reformulated 
higher-order theory results match the exact analytical results very closely, both qualitatively and 
quantitatively. The normal stress contour plots obtained from the reformulated higher-order the- 
ory are perfectly symmetric about the horizontal and the vertical lines passing through the center of 
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Fig. 16. Boundary conditions and mesh discretization for the Eshelby problem. 
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Fig. 17. Comparison of stress fields obtained using HOTFGM and analytical solutions. 
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the fiber inclusion as they ought to be due to the symmetry of the problem about these lines. Also, 
the normal stresses are nearly uniform in the fiber and the shear stress is nearly zero. The stress 
fields in the vicinity of the fiber are almost the same as those of the exact analytical solution. 

The higher-order theory is sometimes confused with the finite-element technique. To demon- 
strate that the two approaches are fundamentally different, a finite-element solution was developed 
for the above problem. The same mesh was generated in ANSYS using 2-D structural 4-node 
(Planel82) and 2-D structural 8-node (Planel83) elements with unit negative pressure applied along 
X 2 = 0, 1. Planel82 elements use a bilinear approximation of the displacement field while Planel83 
elements have the same order (quadratic) displacement field approximation as in the reformulated 
higher-order theory. The middle nodes along x 2 = 0, 1 and x% = 0, 1 were constrained from moving 
in the x% and X 2 directions, respectively, in order to prevent rigid body motion. Stress contour plots 
generated using the reformulated higher-order theory and finite-element analysis based on the 8- 
noded elements are compared in Fig. 18. The results are plotted using the same color scale and are 
magnified by 250 percent in order to compare the stresses in the vicinity of the fiber inclusion more 
closely. As observed in these plots, the finite-element solution picks up local stress concentrations 
at the interfacial subcell corners while, as seen in Fig. 17, the reformulated higher-order theory 
does not pick these stress concentrations and compares well with the actual solution. Further, the 
traction quantities are not continuous along the element interfaces in the finite-element case while 
they are continuous along the subcell interfaces in the reformulated higher-order theory. Finally, 
the stress components are not fully uniform within the fiber in the finite-element case. Table 2 lists 
the maximum and the minimum stresses obtained from the analytical, reformulated higher-order 
theory and the finite-element solutions. As observed in Table 2, the stress concentrations picked 
in the finite-element analysis based on the 4-noded elements are smaller than the ones picked up 
by the 8-noded elements and hence closer to the actual solution. However, the 4 -noded element 
assumes a bilinear variation of displacement held while the 8-noded element assumes a quadratic 
variation. Hence, the results from the 8-noded element based analysis should be compared with 
the actual solution because of the higher-order held approximation. In the hnite-element case, 
these stress concentrations are picked up because the circular inclusion has been approximated by 
a stair-case pattern shown in Fig. 16 (b). 


(MPa) 

Analytical 

HOTFGM 

FEA 

(4-node) 

FEA 

(8-node) 

-max 

°22 

1.49 

1.58 

2.52 

2.94 

min 

°22 

0.05 

-0.027 

0.02 

0.08 

-max 

a 33 

0.68 

0.72 

0.68 

0.64 

min 

a 33 

-0.35 

-0.3 

-0.63 

-1.65 

-max 

°2Z 

0.37 

0.36 

0.25 

0.72 


-0.37 

-0.36 

-0.25 

-0.72 


Table 2. Comparison of maximum and minimum stresses obtained using analytical, HOTFGM 

and FEA approaches. 

Further, in the hnite-element case, the continuity of displacements is satished in a point-wise 
manner (at the nodes) and therefore, the exact details of the geometry are important in order to 
obtain converged results. This was demonstrated by generating the exact circular inclusion shape 
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Fig. 18. Comparison of stress fields obtained using HOTFGM and FEA. 
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Fig. 19. Comparison of stress fields obtained using FEA with exact circular inclusion and 

analytical solutions. 


NASA/CR— 2002-21 1909 


55 





using 8-noded (Planel83) elements in ANSYS. A total of 1045 elements was used in the analy- 
sis. The results obtained for this case are plotted in Fig. 19 and compared with the analytical 
results. As seen in Fig. 19, the results obtained using the finite-element analysis for the exact 
inclusion geometry match the actual analytical results both qualitatively and quantitatively. In the 
reformulated higher-order theory, the displacement/traction continuity conditions are applied in a 
surface-averaged sense and therefore the various field quantities are blurred acrss the interfaces. 
An approximation of the inclusion geometry produces good results. This is very useful in analyzing 
practical problems where modeling the exact geometry can be very demanding. 

4.4 Validation: Combined Thermomechanical Case 

In this subsection, a combined thermomechanical case is considered. The effect of temperature 
increase on a constrained geometry is studied and the results obtained from the reformulated 
higher-order theory are compared with finite-element analysis. The investigated geometry and the 
thermal boundary conditions are the same as those considered in Subsection 4.1 and are shown 
again in Fig. 20 (a) for convenience. The block is considered to be made up of aluminium and the 
material properties used for the analysis are listed in Table 3. The mechanical constraints imposed 
are the fixidity of the left and the right face of the block. 


Material 

k ( W/m - °C) 

E ( GPa ) 

V 

a x 10“ e (/°C) 

Aluminum 

220 

72.4 

0.33 

22.5 


Table 3. Material properties used for the combined thermomechanical analysis. 

The cross-section of the block in the x 2 — £3 plane was initially discretized into 32 x 32 subcells 
of uniform size as shown in Fig. 20 (b). This discretization was shown to be satisfactory in dealing 
with the thermal problem in Subsection 4.1. As the gradients of various mechanical field quantities 
are expected to be higher around the corners, the mesh should be more refined in the regions of 
higher gradients. Therefore, another discretization with the same number of subcells but with the 
mesh refined around the corners was also considered as shown in Fig. 20 (c). In order to compare 
with finite-element results, the same mesh (32 x 32, non-uniform, Fig. 20 (d)) was created in 
ANSYS. 

The temperature, displacement and stress field contour plots obtained using the reformulated 
higher-order theory and finite-element analysis for the graded (refined around the corners) mesh 
are shown in Figs. 21 and 22. As observed from the contour plots, the temperature and the 
displacement field obtained using the two approaches are visually indistinguishable. The stress 
field contour plots generated using the two methods also match very closely except at the corners. 
The maximum and minimum stresses are obtained at the corners and their magnitudes predicted 
by the reformulated higher-order theory and the finite-element analysis are different as shown in 
Table 4. No consistent trend is observed from these magnitudes. The normal stresses predicted 
by finite-element analysis are higher than the higher-order theory results, while the shear stress 
predicted by the higher-order theory is higher than the finite-element result. A thorough ANSYS 
analysis showed that these stress magnitudes (at the corners) do not converge with any amount of 
mesh refinement but continue to increase as the mesh is refined further and further. This is most 
likely due to the singular behavior at the corners which occurs because of the applied mechanical 
boundary conditions. 
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Problem Definition 



32 x 32 subcells (HOTFGM) 32 x 32 subcells(FEA) 



0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 


(C) (d) 

Fig. 20. Problem definition, investigated geometry and discretization used for validation studies: 

thermomechanical problem. 
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Temperature (HOTFGM) 



0 0.2 0.4 0.6 0.8 



*2 


*2 



(e) 



(f) 


Fig. 21. Temperature and displacement field comparison obtained using HOTFGM and FEA for 
the homogeneous plate with fixed boundary conditions. 
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*3 


(C) 

cr 23 (HOTFGM) 



*2 


FEA 



(b) 


FEA 



(f) 


Fig. 22. Stress field comparison obtained using HOTFGM and FEA for the homogeneous plate 

with fixed boundary conditions. 
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(MPa) 

HOTFGM 

FEA 

,max 

a 22 

-83 

-70 

min 

a 22 

-1871 

-2063 

,max 

ct 33 

152 

59 

min 

ct 33 

-715 

-1253 

,max 

ct 23 

963 

608 

^ n 

-963 

-608 


Table 4. Comparison of maximum and minimum stresses obtained using HOTFGM and FEA 

approaches. 


The various field quantities (temperature, displacements 112 and ^3, and stresses (722,^33, and 
(J23) were also plotted along the cross-section x% = 0.1 in order to demonstrate the advantages of 
refining the mesh in the regions of high gradients. The effect of uniform and non-uniform mesh 
is studied and the results are compared with the finite-element results in Fig. 23. From these 
cross-sectional plots, it is observed that the temperature and displacements converge with 32 x 32 
subcells even for the uniform mesh case. However, the stresses are different for the uniform mesh 
case, especially near the corners. This is because of the higher gradients in these regions. The 
stresses are related to the derivatives of the displacement field and, therefore, the error in the stress 
field is expected to be higher than in the displacement field. The non-uniform mesh with 32 x 32 
subcells converges better then the uniform 32 x 32 subcell mesh. This shows the advantages of 
effectively graded meshes. Therefore, an effective mesh should be created with relatively more 
subcells in the regions of high field gradients. As observed from the above cross-sectional plots, 
comparable results are obtained from the finite-elements and the reformulated higher-order theory 
for the graded mesh. 

5 Application: Thermal Barrier Coatings 

As described in Section 1, one of the most important applications of functionally graded materials is 
in the thermal protection systems such as thermal barrier coatings. In many practical applications 
such as aerospace engines, electronic circuit boards, packaging of chips, etc., the structures are 
subjected to very high thermal gradient loading. The metallic part in these structures yields when 
subjected to very high temperatures. In order to prevent this yielding, the metallic substrate is 
coated with a low conductivity ceramic layer to reduce the temperature to which the metal is 
exposed. These ceramic coatings of metallic substrates with an adhesive layer between the metal 
and ceramic regions, or the bond coat, are called thermal barrier coatings or TBCs. 

Layers of bond coat and ceramic with sharp interfaces produce high interlaminar stresses at the 
edges which may lead to separation of these layers. Therefore, in order to prevent these free-edge 
interlaminar stresses, which occur due to large material property mismatch, the microstructure is 
gradually varied. There are various spray techniques such as plasma spray, flame spray, arc spray, 
etc., which are used in different industries for fabricating thermal barrier coatings. The metallic 
substrate is first sprayed with an adhesive which acts as a bond coat between the metal and the 
ceramic material, followed by the coating itself. These spray techniques can be very efficient in 
producing functionally graded microstructures for thermal barrier coatings. 

Thermal barrier coatings can be analyzed by using either the uncoupled or the coupled approach 
as described in Section 1. In the uncoupled approach, the gradually varying microstructure is either 
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Temperature field T 


Stress field a 



Fig. 23. Comparison of temperature displacement and stress distributions along X 3 = 0. 

obtained using HOTFGM and FEA. 
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assumed to be continuous function of spatial coordinates which best fits the given distribution or 
approximated by layers with constant homogenized properties. These homogenized properties are 
calculated using various micromechanical models. However, the uncoupled approach neglects the 
locally produced effects of microstructural gradation. Using the layered approach with homogenized 
properties may produce stress concentrations at the interfaces in the vicinity of free edge which 
would not otherwise occur in the actual microstructure. 

In this section, a TBC with continuously varying microstructure is analyzed using the reformu- 
lated higher-order theory and the results are compared with finite-element analysis. 

5.1 Problem Definition 

The investigated geometry, microstructural gradation and the boundary conditions are shown in 
Fig. 24 (a). The cross-section of the TBC in the X 2 — X 3 plane has unit dimension in the x 3 direction 
and is twice as long in the X 2 direction. The cross-section is discretized into 120 x 60 subcells such 
that each subcell has an aspect ratio of one. The dimension of the block is considered infinite in 
the out of plane (x\) direction and plane strain analysis is carried out. The thermal barrier coating 
is subjected to zero temperature at the bottom face and a concentrated temperature at the top. 
The left and the right faces are insulated against conduction (zero heat flux). At the top face, the 
temperature is assumed to vary with the x 2 coordinate according to the exponential function 

T(x 2 , 1) = 1325 (cos \x 2 ~ 1|) 20 + 25 (207) 

as shown in Fig. 24 (a). The temperature is calculated at the right corner of the top edge of each 
subcell and that constant value is applied at the top face of the subcell. The bottom face of the TBC 
is placed on rollers and the middle two subcells at the bottom are fixed in order to prevent rigid body 
motion. The bottom 12 rows of subcells (/? = 1, 2, .., 12, 7 = 1,2,..., 120) are assigned the properties 
of steel which is the substrate. The next 8 rows of subcells (/? = 13, 14, .., 20, y = l,2, ..., 120) are 
assigned the properties of the alloy CoCrAlY which acts as the bond coat adhesive. In the remaining 
40 rows of subcells (/? = 21,22, ..,60, 7 = 1, 2, ..., 120), the material properties are assigned such 
that the adhesive (CoCrAlY) lies at the bottom and the ceramic (zirconia) at the top with gradual 
variation from the bond coat to the ceramic top coat. The material properties used in the analysis 
are listed in Table 5. The volume fraction of the CoCrAlY bond coat is plotted as a function of x% 
in Fig 24 (b), showing the gradually changing microstructure. 


Material 

k ( W/m - °C) 

E ( GPa ) 

V 

a x 10“ 6 (/°C) 

Steel 

60.5 

207 

0.33 

15 

CoCrAlY 

2.42 

197 

0.25 

11 

Zirconia 

0.5 

36 

0.2 

8 


Table 5. Material properties used for the thermal barrier coating application. 


5.2 HOTFGM and Finite-Element Comparison 

For the given geometry, microstructure and boundary conditions, the various field quantities (tem- 
perature, displacements U 2 and ^3, stresses 022, ^33 and a 2 s) were generated using the reformulated 
higher-order theory and are shown in Figs. 25 (a) - 30 (a). The various field quantities were also cal- 
culated using the finite-element method in order to compare the reformulated higher-order theory 
results. The finite-element solution was obtained by generating the same mesh and microstructure 


NASA/CR— 2002-21 1909 


62 





<-T = 25 




(b) 


Fig. 24. Problem defintion, investigated geometry and material gradation used for TBC 

application. 
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in ANSYS using 2-D coupled field, 4-node (Planel3) elements. In the finite-element case, the 
temperature applied at the top face is calculated at the nodes using the assumed exponential 
function. The bottom face is constrained to move in the X 3 direction and the middle node at the 
bottom is fixed in order to prevent rigid body motion. The generated contour plots for the various 
field quantities are shown in Figs. 25 (b) - 30 (b). Also, the maximum and the minimum values of 
the various field quantities are listed in Table 6. 



Maximum 

Minimum 


HOTFGM 

FEA 

HOTFGM 

FEA 

u 2 

0.0023 

0.0022 

-0.0023 

-0.0022 

113 

0.0035 

0.0034 

0 

0 

022 

103.7 

107.9 

-140.8 

-225.6 

033 

71.5 

105.4 

-186.6 

-217.2 

023 

70.1 

86.7 

-70.9 

-84.7 


Table 6. Comparison of maximum and minimum displacements and stresses obtained using 

HOTFGM and FEA approaches. 

As observed in the contour plots, Figs. 25 - 27, and also in Table 6, the temperature and dis- 
placement fields obtained from the reformulated higher-order theory and the finite-element analysis 
match very closely. However, the stress values obtained using the two approaches are considerably 
different in the graded region, eventhough the general distributions are similar. In the finite- 
element case, the stresses across the interfaces of the elements vary considerably as observed from 
the contour plots in Figs. 28 (b) - 30 (b). The a 22 and a 2 3 stress components are the traction 
components along the vertical interfaces of the subcells and a 23 and (J33 are the traction compo- 
nents along the horizontal interfaces of the subcells. These traction quantities are expected to be 
approximately continuous from one subcell or element to another. In the reformulated higher-order 
theory, the traction quantities are continuous across the subcell interfaces. In the finite-element 
case, on the other hand, the difference in these tractions across the interfaces of the elements, 
particularly in the graded region, is very high, which demonstrates that the stress field has not con- 
verged. This is partly due to the continuity of tractions/displacements which is explicitly applied 
(in a surface- average sense) across the subcell interfaces in the reformulated higher-order theory. 
In the finite-element case, continuity of displacements is explicitly applied in a pointwise manner 
(at the nodes) and then the potential energy is minimized in order to obtain the displacement and 
stress fields. Therefore, in the finite-element case, a more refined mesh is required in order to get 
converged solution for the stresses. Also, as observed in the stress contour plots, the boundary 
conditions are somewhat better satisfied by the reformulated higher-order theory. 

6 Conclusions and Future Work 

6.1 Conclusions 

An efficient reformulation of the higher-order theory for functionally graded materials for two- 
dimensional thermoelastic problems has been successfully developed. The subcell microvariables, 
which were the basic unknowns in the original higher-order theory, were expressed in terms of the in- 
terfacial surface-averaged quantities (temperatures and displacements) and these interfacial surface- 
averaged quantities were considered to be the basic unknown quantities in the reformulation. The 
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X 3 Temperature 



(a) HOTFGM 



(b) FEA 



0 200 400 600 800 1000 1200 


Fig. 25. Temperature field for TBC application obtained using HOTFGM and FEA. 
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(b) FEA 
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Fig. 26 . Displacement field 112 for TBC application obtained using HOTFGM and FEA. 
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(a) HOTFGM 



(b) FEA 



0 0.5 1 1.5 2 2.5 3 


Fig. 27. Displacement field 113 for TBC application obtained using HOTFGM and FEA. 
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(a) HOTFGM 



(b) FEA 



Fig. 28. Stress field <722 for TBC application obtained using HOTFGM and FEA. 
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(a) HOTFGM 



(b) FEA 



Fig. 29. Stress field (733 for TBC application obtained using HOTFGM and FEA. 
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(a) HOTFGM 



(b) FEA 



Fig. 30. Stress field a 23 for TBC application obtained using HOTFGM and FEA. 
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local-global conductivity and local-global stiffness matrix approach was used in order to eliminate 
redundant equations. This resulted in the reduction of the size of the global conductivity and 
stiffness matrices by approximately sixty percent. The reduction in the number of equations has 
enhanced the theory’s capability to analyze computationally intensive and demanding cases. Prob- 
lems requiring greater mesh discretization (large number of subcells) which could not be analyzed 
using the original higher-order theory because of the large number of equations involved are now 
solvable with the reformulated higher-order theory. 

Mesh sensitivity and validation studies were carried out for various thermal, mechanical and 
combined thermomechanical cases. A practical application of the thermal barrier coating was also 
analyzed. The results were compared to analytical and finite-element solutions for the various cases. 

From the mesh sensitivity studies, it was observed that the temperature field converges very 
quickly with mesh refinement and the temperature distribution can be accurately captured using 
a relatively coarse mesh, especially for problems involving homogeneous materials. Also, averaging 
the temperature across subcell interfacial lines calculated using subcell microvariables on each side 
of the interface produces better results. 

A functionally graded case with thermal conductivity varying according to an exponential func- 
tion was considered and thermal analysis was carried out. The results obtained from the higher- 
order theory matched the finite-element analysis very closely. However, in the finite-element case 
the temperature field did not converge to the actual temperature field at the points of tempera- 
ture discontinuity. A more refined mesh is required at these points in the finite-element case. On 
the other hand, the results obtained from the reformulated higher-order theory exhibited better 
convergence at these points. 

In many practical problems, the mesh discretization required to simulate exact inclusion shapes 
can be very demanding. The Eshelby problem demonstrated the efficiency of higher-order theory 
in terms of mesh discretization approximation of the inclusion phase. It was observed that a 
decent approximation of the circular inclusion using a rectangular grid produced results which were 
comparable with the actual analytical solution. However, in the finite-element case, approximating 
an inclusion shape using the same mesh discretization picks up stress concentrations at the sharp 
edges which are not present in the actual solution. Therefore, a more refined mesh is required in 
order to properly approximate the shape of inclusions in the finite-element case. 

A thermomechanical case involving a homogeneous plate was also analyzed. The results from 
the higher-order theory and the finite-elements matched very closely. The efficiency of using a 
locally refined mesh in the regions of high temperature and stress gradients was demonstrated. 

In the thermal barrier coating application, it was observed that the temperature and displace- 
ment fields obtained using the reformulated higher-order theory and the finite-element approach 
were visually indistinguishable, in contrast to the stress fields. In the finite-element case, the 
tractions were not continuous across the element interfaces, especially in the graded region, and 
the boundary conditions were not as well satisfied. In the reformulated higher-order theory, on 
the other hand, the tractions were approximately continuous across the subcell interfaces and the 
boundary conditions were better satisfied. This demonstrated the efficiency of the higher-order 
theory in analyzing functionally graded microstructures. In the case of functionally graded mate- 
rials analyzed using the finite-element approach, mesh discretization required for converged results 
can be computationally very demanding if the actual microstructural details are explicitly taken 
into account. In the finite-element case, displacement continuity is satisfied in a point wise man- 
ner and hence stress concentrations are picked up at the points of material discontinuity at sharp 
edges. However, in the higher-order theory, the continuity of tractions/displacements is applied 
in a surface-average sense. This smoothing operation produces sufficiently accurate solutions with 
relatively coarse meshes and approximate inclusion shapes. 
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6.2 Future work 


The completed two-dimensional thermoelastic reformulation of the higher-order theory is useful for 
plane strain and generalized plane strain problems. In many practical applications, the type of 
problems encountered are not limited to two dimensions, however. Many cases involve out-of-plane 
loading conditions. Therefore, the next step is to extend the reformulation to three dimensions in 
order to analyze more general problems. Also, in cases involving traditional composites, the material 
is locally heterogeneous but globally periodic. Structural problems involving composites can be 
analyzed by homogenizing the material using the concept of a representative volume element or 
RVE. The effective properties of the locally heterogeneous RVE can be found using the homogenized 
version of the reformulated higher-order theory by applying periodic boundary conditions, and 
then the structure can be analyzed globally using the reformulated higher-order theory or other 
techniques. The second step, therefore, is to incorporate periodic boundary conditions and out-of- 
plane loading capabilities. 

Generally, in practical problems encountered in the industry, deformations enter into the inelas- 
tic range at the point where the structure can still withstand considerably higher loads. Also, the 
loading history affects stresses and deformations to an extent that depends on the type of material. 
In order to carry out an optimum design of a structure, it is necessary to include these plastic 
and viscoelastic or viscoplastic effects and analyze the design taking into account such behavior. 
Therefore, viscoelastic, viscoplastic and plastic capabilities should be included in the reformulated 
higher-order theory. 

Presently, material properties are considered to be constant within a subcell. The material 
properties can be assumed to be linear or higher-order in local coordinates depending on the 
actual gradation of the microstructure. This may potentially lead to further reductions in mesh 
discretization required to obtain converged solutions. The capability of having linearly varying 
material properties within a subcell has been incorporated by Zhong (2002) into the thermal portion 
of the reformulated higher-order theory. This will be extended to the mechanical portion. 

The stringent environment and loading conditions to which structures are subjected may lead 
to the initiation of cracks. Another long-term goal is to include fracture mechanics capabilities in 
the reformulated higher-order theory in order to be able to analyze crack problems. 
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7 Appendix 
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The submatrix used in the global stiffness matrix assembly 
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